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Abstract 

This  thesis  proposes  and  analyzes  a  new  measure  of  functional  capability  for 
satellite  constellations  that  incorporates  the  instantaneous  availability  and  mission 
effectiveness  of  individual  satellites.  The  capability  measure  yields  a  continuous  score 
between  zero  and  one  accounting  for  the  degree  to  which  the  constellation  meets 
operational  requirements.  The  measure  is  computed  from  an  average  of  satellite  ca¬ 
pabilities,  composed  of  the  product  of  the  satellite’s  instantaneous  availability  and 
value  score.  Instantaneous  availability  is  acquired  by  modeling  the  satellite  degrada¬ 
tion  status  as  either  a  time-homogenous,  continuous-time  Markov  chain  (CTMC)  if  it 
possesses  functions  with  exponential  lifetime  distributions,  or  as  a  time-homogenous, 
semi-Markov  process  (SMP)  if  the  function  lifetime  distributions  are  not  exponen¬ 
tial.  The  satellite  value  score  represents  the  individual  satellite’s  contribution  to  the 
overall  constellation  mission  and  is  obtained  using  multi-attribute  value  theory.  For 
illustrative  purposes,  analytical  results  were  compared  with  those  obtained  via  the 
Monte  Carlo  method  and  were  found  to  be  indistinguishable  with  substantially  less 
computational  effort. 
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STOCHASTIC  CAPABILITY  MODELS 
FOR  DEGRADING  SATELLITE  CONSTELLATIONS 


1.  Introduction 

1.1  Background 

The  unique  vantage  point  of  space  offers  numerous  advantages  to  military  forces 
operating  in  a  global  arena.  Missions  employed  by  the  United  States  military  in  the 
present,  and  in  the  future,  will  rely  heavily  on  information  provided  by  space  based 
systems.  A  minimum  level  of  capability  is  required  to  support  these  missions,  and 
the  efficient  management  of  maintaining  this  capability  level  into  the  future  depends 
on  an  accurate  capability  assessment  via  performance  and  availability  measures. 
Without  an  idea  of  how  well  these  systems  meet  the  requirements  of  the  missions  that 
rely  on  their  services,  the  potential  exists  for  excessive  spending  to  ensure  systems  are 
capable,  or  systems  may  not  be  adequately  capable  to  support  military  operations. 
Justification  for  funding  new  systems  and  maintaining  current  systems  may  occur 
due  to  the  gap  that  exists  between  current  capabilities  and  the  requirements  of  the 
user. 

Missions  performed  by  U.S.  military  satellite  constellations  are  summarized 
in  the  Joint  Doctrine  for  Space  Operations  [5]  published  in  August  2002.  These 
satellite  constellations  perform  many  essential  missions  pertaining  to  intelligence, 
surveillance,  and  reconnaissance  (ISR);  integrated  tactical  warning  and  attack  as¬ 
sessment;  environmental  monitoring;  communications;  and  position,  velocity,  time, 
and  navigation.  The  increased  importance  of  these  missions  to  the  successful  com¬ 
pletion  of  military  goals  and  the  increased  likelihood  of  attacks  on  space  assets  which 
support  these  missions  has  generated  a  need  for  an  accurate  measurement  of  a  satel- 
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lite  constellation’s  functional  capability,  a  measure  of  the  system’s  ability  to  perform 
its  intended  mission. 

Deriving  a  functional  capability  methodology  for  satellite  constellations  is  of 
interest  to  any  company  or  agency  that  provides  direct  support  to  satellite  constel¬ 
lations.  Capability-based  constellation  replenishment  models  require  an  accurate 
assessment  of  a  constellation’s  capability.  This  functional  capability  methodology 
may  also  help  redefine  maintenance  and  operation  policies,  increasing  the  life  ex¬ 
pectancy  of  the  constellation  by  perhaps  operating  its  satellites  in  a  more  efficient 
manner  relative  to  the  satellite’s  current  mission  and  condition.  Also,  any  sudden 
change  in  availability  of  the  constellation  may  indicate  if  an  attack  on  space  assets 
has  occurred  and  if  so,  when  the  attack  occurred  and  to  what  degree  the  attack  had 
on  the  system. 

The  incorporation  of  new  technology  and  the  sustainment  of  current  capabili¬ 
ties  involving  space  systems  is  paramount  to  the  success  of  future  military  operations. 
Joint  Vision  2020  [6]  is  a  collaborative  strategy  composed  by  the  Joint  Chiefs  of  Staff 
that  describes  the  necessary  transformations  of  the  joint  military  forces  current  ca¬ 
pabilities  for  successful  future  military  and  humanitarian  operations.  The  document 
states  two  important  philosophies  relative  to  military  operations  in  space.  First, 
information  superiority,  “the  capability  to  collect,  process,  and  disseminate  an  un¬ 
interrupted  flow  of  information  while  exploiting  or  denying  an  adversary’s  ability 
to  do  the  same  [6:8],”  is  recognized  as  a  necessary  component  in  successful  military 
operations.  The  publication  emphasizes  the  adoption  of  a  new  doctrine  dependent 
on  the  availability  of  timely  and  accurate  information-  information  assumed  to  be 
readily  available  in  a  real  time,  error-free  delivery.  Second,  in  order  to  be  successful 
against  enemies  of  the  future,  military  capabilities  will  have  to  be  adequately  dy¬ 
namic  to  keep  the  enemy  from  being  able  to  adapt  and  counteract  our  capabilities. 
The  vision  statement  suggests  that  one  of  the  best  ways  to  keep  military  capabilities 
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from  entering  a  static  state  is  through  “...the  steady  infusion  of  new  technology  and 
modernization  and  replacement  of  equipment”  [6:3]. 

1.2  Problem  Definition  and  Methodology 

The  United  States  Air  Force  has  a  requirement  to  measure  satellite  perfor¬ 
mance  based  on  the  capability  of  a  satellite  constellation,  rather  than  measures 
derived  from  expected  lifetime  estimates  [9] .  This  thesis  will  propose  a  model  which 
describes  satellite  functional  capability  on  a  continuous  scale  and  can  assess  the 
partial  capability  of  a  constellation  assuming  certain  conditions. 

As  with  all  environments  in  which  the  military  operates,  space  offers  unique 
advantages  and  disadvantages  that  affect  system  functional  capability.  The  obvious 
advantages  to  military  systems  include  global  access,  the  ability  to  have  a  line  of 
sight  to  any  point  on  the  globe,  and  persistence ,  the  ability  of  space  assets  to  stay  on 
orbit  for  long  durations  of  time.  There  are  also  many  disadvantages  to  operating  in 
space  that  can  have  an  effect  on  functional  capability  of  these  systems.  First,  space  is 
a  hostile  environment  which  can  cause  system  failures,  disruptions,  and  degradation 
in  performance,  all  of  which  directly  affect  the  system’s  state  of  availability.  Second, 
the  location  of  space  assets  limits  convenient  access  to  these  systems,  resulting  in 
high  costs  associated  with  fielding  and  maintaining  them.  Also,  the  amount  of 
expendable  resources  these  systems  are  designed  to  store  is  restricted  due  primarily  to 
deployment  restrictions  associated  with  weight  and  the  extent  to  which  maintenance 
can  be  performed  to  replenish  degraded  systems. 

Historically,  metrics  based  solely  on  lifetime  estimates  have  been  used  in  place 
of  capability-based  measures  for  satellites,  which  were  then  pooled  into  a  measure 
representing  the  constellation  lifetime.  Mean  Mission  Duration  (MMD),  the  expected 
lifetime  of  the  space  asset,  was  used  to  estimate  how  long  the  system  would  be 
available.  This  measure,  however,  does  not  inherently  incorporate  how  well  the 
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system  performs  its  mission.  In  most  cases,  MMD  underestimates  the  lifetime  of  the 
system  and  lacks  the  ability  to  accurately  predict  the  capability  of  the  system  at  a 
future  point  in  time  [9].  Satellite  systems  are  over-engineered  to  meet  guaranteed 
satellite  lifetime  requirements.  Derating  is  a  method  used  to  increase  the  reliability  of 
a  component  by  operating  it  at  a  significantly  lower  stress  level  than  it  was  originally 
designed.  Higher-rated  components  are  used  in  place  of  normal  rated  components 
to  increase  the  reliability  and  lifetime  of  the  component.  Over-engineering  leads  to 
systems  that  meet  their  required  MMD  specification  with  a  high  probability,  making 
MMD  a  poor  measure  upon  which  to  base  an  estimate  of  the  system’s  true  lifetime 
[9], 

Most  current  methods  of  assessing  constellation  capability  use  either  decision 
analysis  techniques  or  reliability  measures  using  steady-state  system  availability  as 
the  basis  of  their  methodology.  Multi- attribute  value  theory,  a  subset  of  decision 
analysis,  may  provide  a  methodology  of  scoring  a  system  by  assigning  weights  to 
certain  attributes  of  the  system  directly  correlated  to  a  defined  value,  in  the  case  of 
functional  capability,  a  value  representing  the  system  performing  its  mission.  The 
shift  from  alternative  based  thinking  to  value  focused  thinking  provides  an  approach 
to  assessing  the  value  of  a  satellite  to  the  overall  constellation  mission  and  may 
help  in  generating  the  overall  measure  of  functional  capability.  The  process  involves 
first  determining  the  appropriate  attributes,  usually  through  interviews  with  subject 
matter  experts,  assigning  a  score  related  to  order  of  priority  of  these  attributes, 
then  determining  a  weight  of  importance  for  each  attribute.  After  attributes  are 
determined  with  their  corresponding  values  and  weights,  a  function  is  composed 
which  yields  a  value  between  0  and  1  that  indicates  how  well  the  system  is  able  to 
perform  its  mission.  This  approach  yields  an  estimate  of  functional  capability  and 
is  strictly  confined  to  a  particular  constellation,  since  the  attribute  selection  and 
scoring  and  weighting  are  dependent  on  the  nature  of  a  constellation  mission. 
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The  availability  of  a  system  is  defined  in  terms  of  its  reliability  and  maintain¬ 
ability,  all  of  which  are  probability  functions  of  time.  Reliability  is  the  probability 
that  a  system  does  not  fail  before  a  specific  time;  it  is  the  complement  to  the  probabil¬ 
ity  that  a  system  fails  before  a  specific  time,  also  known  as  the  cumulative  probability 
distribution  function.  Maintainability  is  the  probability  that  a  system’s  reliability 
can  be  restored  to  some  level  through  maintenance  or  repair  actions.  Instantaneous 
availability  is  the  probability  a  system  is  available  to  perform  its  designed  mission 
at  some  point  in  time.  Steady-state  availability  is  the  proportion  of  up  time  to  the 
total  time  of  operation  over  a  very  long  period  of  time.  If  a  system’s  steady-state 
availability  exists,  the  system’s  instantaneous  availability  will  eventually  converge  to 
its  steady-state  availability. 

System  effectiveness  is  a  measure  which  gives  the  probability  that  a  system 
will  perform  its  intended  mission  under  its  designed  operating  conditions.  Ebel- 
ing  [14],  suggests  that  system  effectiveness  is  comprised  of  the  system’s  operational 
readiness,  the  probability  of  the  system  working  at  the  start  of  its  life;  the  system’s 
mission  availability,  the  proportion  of  lifetime  the  system  is  performing  its  mission; 
and  the  system’s  design  adequacy,  the  probability  the  system  will  perform  as  it  was 
designed.  In  terms  of  satellite  effectiveness  in  a  constellation,  the  operational  readi¬ 
ness  is  described  as  the  probability  the  satellite  is  on  station  and  working,  surviving 
the  launch,  deployment,  and  initialization  processes;  mission  availability  is  described 
as  the  proportion  of  the  satellite’s  lifetime  in  the  constellation  performing  its  mission; 
and  the  satellite’s  design  adequacy  is  described  as  the  probability  the  satellite  does 
what  it  was  designed  to  do.  If  operational  readiness,  mission  availability,  and  design 
adequacy  are  independent  probabilities,  one  possible  definition  of  satellite  effective¬ 
ness  might  be  the  product  of  its  respective  probabilities  of  operational  readiness, 
mission  availability,  and  design  adequacy. 

The  effectiveness  measure  results  in  a  number  between  0  and  1  in  a  multi¬ 
plicative  manner,  which  describes  a  level  of  effectiveness  for  the  constellation  with 
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no  associated  units.  The  system  effectiveness  measure  is  based  on  the  engineering 
side  of  system  development,  where  design  adequacy  is  a  percentage  based  on  speci¬ 
fications  that  are  more  focused  on  design  versus  operational  adequacy.  Ebeling  [14] 
suggests  that  one  might  obtain  the  design  adequacy  measure  by  finding  the  propor¬ 
tion  of  time  the  system  adequately  performs  its  intended  job.  However,  a  subject 
matter  expert  would  need  to  be  involved  to  make  the  distinction  between  what  cri¬ 
teria  to  use  in  measuring  adequacy.  Decision  analysis  techniques  could  then  be  used 
to  assess  how  well  the  system  performed  its  job. 

Efforts  have  been  made  to  create  a  meaningful  estimate  of  constellation  func¬ 
tional  capability  which  contain  elements  of  reliability  techniques  combined  with  deci¬ 
sion  analysis  techniques  [9].  A  score  is  formed  by  taking  the  product  of  three  terms: 
the  product  of  the  reliability  associated  with  the  random  and  wear-out  phases  of 
the  satellite’s  lifetime,  which  accounts  for  the  probability  of  the  satellite  being  on 
station  and  currently  working  at  a  specific  time;  the  duty  cycle  ratio,  the  ratio  of 
estimated  mean  duty  cycle  over  the  beginning-of-life  duty  cycle,  which  addresses 
the  percent  of  time  the  satellite  is  performing  its  mission;  and  the  band  capabil¬ 
ity,  which  is  a  formulation  of  a  subject  matter  expert’s  assessment  of  the  satellite’s 
payload  capability. 

A  satellite’s  duty  cycle  is  the  proportion  of  time  the  satellite  is  performing  its 
mission.  When  a  satellite  is  in  orbit,  there  may  be  cycles  of  time  in  which  there 
are  no  tasks  to  be  performed  by  the  satellite.  Changing  the  mode  of  the  satellite 
to  a  standby  status  during  periods  of  inactivity  extends  the  overall  lifetime  of  the 
satellite.  Duty  cycle  is  directly  related  to  the  reliability  of  the  satellite  over  time; 
the  higher  the  duty  cycle,  the  faster  the  functions  of  the  satellite  will  wear  out. 

The  composite  score  was  designed  in  the  hopes  of  capturing  the  constellation’s 
satellite  equivalence,  an  estimate  of  the  equivalence  of  how  many  fully  capable  satel¬ 
lites  are  in  the  constellation.  It  generates  a  number  which  represents  the  number  of 
satellite  equivalents  in  the  constellation  in  an  additive  manner  where  equivalence  is 
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relative  to  the  capability  level  of  the  newest  satellite.  The  number  is  then  compared 
to  a  minimum  number  of  satellites  required  to  accomplish  a  mission.  Theoretically, 
one  could  define  intervals  of  satellite  numbers  which  describe  states  of  the  constel¬ 
lation’s  capability  to  include  fully  capable,  partially  degraded,  and  fully  degraded 
states. 

The  reliability  measures  for  the  random  and  wear-out  phases  of  the  satellites 
are  formed  by  statistical  models  based  on  historical  observation  data  from  similar 
satellites.  These  models  are  open  to  errors  associated  with  small,  dependent  samples 
from  systems  which  may  not  possess  reliability  characteristics  of  the  systems  that 
are  being  modeled. 

Basing  the  composite  score  on  a  theoretical  number  of  satellites  in  a  constel¬ 
lation  would  be  acceptable  if  the  number  of  satellites  is  a  representative  estimate 
for  the  constellation’s  capability.  This  may  or  may  not  be  a  good  measure  for  some 
constellations  due  to  the  nature  of  their  missions  as  well  as  the  nature  of  their  config¬ 
urations.  For  example,  in  a  GPS  constellation,  numeracy  has  a  definite  correlation  to 
area  coverage  due  to  line-of-sight  restrictions  in  the  transfer  of  data  between  ground 
and  space.  Considering  a  measure  which  represents  a  satellite  equivalency  makes 
sense  for  such  a  constellation  since  the  number  of  satellites  holds  direct  relation  with 
the  constellation’s  capability.  However,  there  may  be  instances  when  this  measure 
may  not  fully  explain  the  satellite  capability.  For  example,  if  there  exist  satellites 
within  the  constellation  that  depend  on  other  relay  satellites  for  communication  with 
ground  stations,  these  relay  satellites  hold  more  value  in  terms  of  the  mission,  and 
a  score  based  on  the  number  of  satellite  equivalents  does  not  differentiate  between 
non-identical  satellites. 

The  goal  of  this  thesis  is  to  present  a  methodology  that  combines  the  capa¬ 
bility  assessment  of  decision  analysis  techniques  with  reliability  theory  to  develop  a 
measure  that  meets  the  requirements  of  the  United  States  Air  Force.  This  measure 
will  reflect  the  ability  of  a  constellation  to  fulfill  mission  requirements  and  it  will 
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encompass  the  probability  associated  with  the  constellation  being  available  to  per¬ 
form  its  mission.  The  methodology  will  assess  the  constellation  at  the  satellite  level, 
where  methods  for  deriving  its  instantaneous  availability  will  be  the  main  focus. 
Using  instantaneous  availability  gives  a  more  representative  measure  of  availabil¬ 
ity  for  newer  constellations  that  have  been  on  orbit  for  a  relatively  short  amount  of 
time.  The  decision  analysis  methodology  will  be  applied  to  the  satellite’s  operational 
contribution  to  the  constellation,  focusing  on  the  extent  to  which  the  satellite  ac¬ 
complishes  constellation  mission  requirements,  rather  than  its  ability  to  meet  design 
specifications. 

1.3  Thesis  Outline 

The  next  chapter  will  provide  an  overview  of  the  literature  which  contributes 
to  the  development  of  functional  capability  measures,  including  basic  concepts  of 
reliability  and  multi-attribute  value  theory  followed  by  a  summary  of  methods  that 
attempt  to  optimize  constellation  replenishment  policies  and  compute  constellation 
performance  measures.  Chapters  3  and  4  will  discuss  the  methodology  behind  the 
new  constellation  functional  capability  measure.  Chapter  5  will  apply  the  methodol¬ 
ogy  to  navigation,  communication,  and  meteorological  satellite  constellation  exam¬ 
ples  to  illustrate  the  methodology  and  to  compare  analytical  results  with  simulation 
results,  demonstrating  advantages  of  the  analytical  methods.  Finally,  Chapter  6  will 
summarize  the  results  and  offer  suggestions  for  further  research. 
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2.  Review  of  the  Literature 

This  chapter  will  discuss  the  literature  that  has  a  direct  contribution  to  the 
understanding  and  development  of  a  satellite  constellation  functional  capability  mea¬ 
sure.  The  first  section  will  present  a  timeline  of  general  developments  in  the  field  of 
reliability  pertaining  to  the  characterization  of  system  lifetimes.  The  second  section 
will  focus  on  the  main  developments  of  availability.  The  third  section  will  discuss 
some  of  the  current  constellation  replenishment  policies  in  the  literature.  The  fourth 
section  will  end  the  chapter  with  a  discussion  on  the  methodologies  related  to  ob¬ 
taining  a  satellite  capability  measure  either  by  applying  the  theory  of  reliability  or 
multi-attribute  utility  theory. 

2. 1  Reliability 

Reliability  theory  is  the  application  of  probability  and  statistics  in  determining 
characteristics  of  a  system’s  failures  over  time.  The  development  of  this  theory 
started  during  our  nation’s  transition  from  an  agricultural  to  an  industrial  based 
economy.  A  need  for  machines  to  be  maintained  in  an  efficient  manner  was  key 
to  a  profitable  manufacturing  process.  One  branch  of  the  development  of  reliability 
focuses  on  the  lifetime  of  systems,  while  another  focuses  on  characterizing  an  optimal 
policy  of  repairing  and  replacing  these  machines.  The  paper  by  Barlow  [2]  gives  a 
history  of  the  main  developments  in  reliability  theory. 

Reliability  is  a  probability  function  of  time  that  indicates  the  likelihood  of  a 
system  not  failing  before  a  certain  time.  This  function  is  sometimes  also  called  the 
survivor  function,  as  it  characterizes  the  probability  the  system  will  survive  up  to  a 
point  in  time.  The  complement  of  the  reliability  function  is  the  lifetime  cumulative 
distribution  function  (c.d.f.). 

The  early  developments  of  reliability  started  in  the  1930s  and  were  based  on 
hireling  reliability  functions  that  described  the  fatigue  life  of  materials.  The  jus- 
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tihcation  of  modeling  incoming  telephone  calls  as  Poisson  arrivals  by  Erlang  and 
Palm  laid  the  foundation  for  arguments  supporting  the  use  of  exponential  functions 
to  describe  system  lifetimes  [2],  In  the  1950s,  Epstein  and  Sobel  [15]  initiated  the 
assumption  that  the  exponential  distribution  was  most  applicable  to  model  system 
reliability,  which  was  later  supported  by  Davis  in  1952  when  he  published  a  paper 
characterizing  the  fit  of  several  reliability  functions  to  actual  failure  data  [12],  The 
exponential  function  was  popular  not  only  because  of  the  research  supporting  it,  but 
also  due  to  the  simple  results  obtained  with  its  use.  In  1939,  Weibull  characterized 
a  different  function  which  described  the  breaking  strengths  of  materials  [40].  This 
reliability  function  was  first  adopted  due  to  its  simplicity  and  was  later  given  more 
attention  by  papers  from  Kao  [27]  and  from  Zclcn  and  Dannemiller  [43]  to  be  more 
robust  than  the  exponential  function  in  certain  applications. 

2.2  Availability 

Availability  is  the  probability  a  system  is  available  to  perform  its  intended 
function  at  some  time  t,  given  that  the  system  may  have  failed  and  received  a  repair 
in  the  past.  It  is  a  measure  that  is  closely  tied  to  the  concept  of  system  reliability 
by  encompassing  the  maintainability  aspect  of  a  system  by  by  accounting  for  the 
increases  in  reliability  caused  by  repairs  or  replacements.  The  availability  measure 
will  always  be  greater  than  or  equal  to  the  reliability  measure;  if  a  system  cannot 
be  repaired  or  replaced,  its  availability  will  be  equal  to  its  reliability.  The  branch  of 
reliability  theory  characterizing  replacement  policies  started  in  the  1940s  with  the 
research  of  Lotka  [31]  and  Campbell  [4], 

The  modeling  of  a  system’s  availability  is  usually  approached  using  a  combina¬ 
torial  model,  a  state-space  model,  or  a  combination  of  the  two.  Combinatorial  models 
include  reliability  block  diagrams,  reliability  graphs,  and  fault  trees.  These  methods 
were  developed  in  the  1970s,  were  the  first  ideas  of  fault  tree  analysis  were  presented 
by  Fussell  and  Vesely  [20]  describing  minimum  cut  sets  for  fault  trees.  The  idea 
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of  combinatorial  models  stems  from  the  assumption  that  the  failure  of  components 
of  a  system  are  independent  of  one  another.  By  arranging  these  components  into 
combinations  of  serial  and  parallel  configurations,  and  assigning  reliability  measures 
to  each  component,  the  system’s  availability  can  be  determined.  By  enumerating  the 
possible  events  that  can  cause  a  system  failure,  the  probability  of  a  system  failure 
can  be  computed  using  the  basic  laws  of  probability. 

State-space  models  allow  the  modeling  of  interactions  between  component  fail¬ 
ures.  A  typical  model  consists  of  describing  the  system  in  terms  of  states  that 
represent  the  condition  in  which  the  system  may  be  found.  The  system’s  degrada¬ 
tion  can  then  be  modeled  as  a  stochastic  process.  The  probability  of  the  system 
being  in  a  state  at  a  given  time  may  be  calculated  based  on  information  on  the  rates 
of  transitions,  provided  the  stochastic  model  possesses  certain  characteristics.  One 
of  the  first  papers  regarding  stochastic  modeling  of  system  reliability  was  written 
by  Weiss  in  1956,  where  he  used  semi-Markov  processes  to  characterize  system  re¬ 
placement  policies  [41].  Two  examples  of  modeling  the  availability  of  systems  using 
state-space  models  are  given  in  papers  by  Fricks  [19]  and  Ibe  [26].  Fricks  discusses 
the  availability  analysis  of  a  multi-computer  system  while  Ibe  describes  the  avail¬ 
ability  modeling  of  a  management  system  which  automates  the  operation  of  a  utility 
company. 

2.3  Satellite  Constellation  Replenishment  Models 

There  are  many  papers  written  on  the  subject  of  satellite  constellation  replen¬ 
ishment.  These  replenishment  models  are  based  on  optimizing  a  certain  aspect  of 
the  constellation,  whether  it  be  minimizing  cost  or  maximizing  capability.  These 
models  depend  on  a  performance  measure  on  which  to  base  the  optimization  model. 

For  example,  in  a  paper  by  Collopy  [8],  a  replacement  model  is  developed  for  a 
satellite  constellation  where  the  problem  of  determining  the  number  of  optimal  spares 
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is  addressed.  The  policy  is  based  on  keeping  the  satellite  constellation  in  a  capable 
state,  defined  by  a  minimum  number  of  operating  satellites.  The  constellation  has 
n  —  m  +  i  satellites  where  m  is  equal  to  the  minimum  number  of  required  satellites, 
and  %  is  equal  to  the  number  of  on-orbit  spares.  The  methodology  starts  by  looking 
at  the  event  in  time  in  which  a  satellite  in  the  constellation  fails.  At  this  point,  a 
spare  replaces  the  failed  satellite  and  the  number  of  satellites  in  the  constellation 
is  reduced  to  n  =  m  +  i  —  1.  The  paper  addresses  the  probability  of  running  out 
of  spares  and  falling  below  the  minimum  number  of  operating  satellites,  n,  before  a 
satellite  replacement  can  be  launched  to  replenish  the  constellation.  This  interval  is 
from  time  of  failure  to  replenishment  of  the  constellation  by  a  non-on-orbit  satellite 
is  referred  to  as  Replacement  Launch.  The  number  of  failures  in  that  interval  has  a 
Poisson  distribution  with  parameter  A  —  n  ■  p,  where  p  is  the  ratio  of  Replacement 
Launch  over  the  MTTF  of  any  given  satellite,  and  n  is  the  number  of  satellites  left 
in  the  constellation;  satellites  are  considered  to  be  independent  and  identical. 

The  probability  of  running  out  of  spare  satellites  before  a  satellite  can  be 
launched  to  replenish  the  constellation  is  equivalent  to  a  system  failure  and  is  given 
as 


which  is  the  probability  of  i  failures,  and  does  not  account  for  any  probability  asso¬ 
ciated  with  more  than  %  failures. 

The  paper  then  defines  the  term  Cycle  as  being  equal  to  satellite  reliability, 
R(t)  divided  by  the  number  of  satellites  in  orbit  prior  to  the  failure,  added  with  the 
time  interval  Replacement  Launch ;  Cycle  is  the  time  interval  in  which  the  probability 
of  system  failure  is  based  on,  where 

R(t) 

Cycle  = - b  Replacement  Launch. 

m  +  i 
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The  system  MTBF  is  given  as  the  Cycle  divided  by  the  probability  of  a  system 
failure  occurring,  where 

MTTFS  =  (  'I'1'  =  Cycle  ■  ex  • 

P(i)  A* 

Collopy  [8]  uses  optimization  techniques  to  minimize  the  number  of  spares 
maintained  in  the  constellation  through  an  objective  function  which  captures  the 
cost  of  the  system  based  on  cost  of  individual  satellites  and  cost  of  a  system  failure. 
The  model  assumes  the  satellite  constellation  is  “capable”  based  on  a  minimum 
number  of  operating  satellites  and  does  not  address  partially  operating  satellites. 
Failed  satellites  are  replaced  by  spares  which  are  either  on  orbit  or  are  waiting  to  be 
launched.  The  satellite  lifetimes  are  modeled  as  exponential  random  variables,  and 
all  satellites  are  independent  and  identical. 

An  example  of  applying  replenishment  methodology  to  satellite  constellations 
is  given  by  Feuchter  [17].  He  discusses  the  conditions  when  a  constellation  should 
be  maintained  by  replacing  failed  satellites  or  by  implementing  on-orbit  repairs  to 
degraded  satellites.  An  optimal  replacement  policy  attempts  to  prescribe  when  it 
is  more  economical  to  repair  the  satellite  versus  simply  replacing  it.  Examples  of 
satellites  designed  for  on-orbit  repair  are  given  as  a  contrast  to  the  types  of  satellites 
the  Department  of  Defense  typically  uses.  For  example,  the  NASA  program  satellites 
such  as  the  Hubble  Space  Telescope  (HST),  the  Gamma  Ray  Observatory  (GRO),  the 
Advanced  X-Ray  Astrophysics  Facility  (AXAF),  and  the  Space  Infrared  Telescope 
Facility  (SIRTF)  are  large,  costly  satellites,  whereas  the  Department  of  Defense 
typically  uses  cheaper,  constellation-based  satellites,  which  are  designed  under  a 
replacement  policy  and  are  deployed  in  orbits  that  have  limited  access. 

Although  the  Department  of  Defense  does  not  plan  on  converting  over  to  a 
repair  policy  in  the  near  future,  Feuchter’s  paper  attempts  to  define  the  boundary 
between  the  two  policies  by  considering  existing  constellations,  culminating  in  the 
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Comprehensive  On-orbit  Maintenance  Assessment  (COMA)  final  report.  “COMA 
focuses  on  constellation  support,  which  consists  of  establishing  a  constellation  and 
keeping  it  functional  throughout  its  life”  [17:499].  It  relies  on  resource  requirements 
generated  by  a  simulation  model,  which  consists  of  a  function  of  constellation  avail¬ 
ability  (differential  lifecycle  cost  model).  Satellite  mass  turns  out  to  be  a  differentia¬ 
tor  between  expendable  and  repairable  satellites  and  as  such,  is  an  important  model 
input. 

The  satellite  is  broken  down  into  six  subsystems;  the  structure,  attitude  control 
system,  communications,  telemetry  tracking  and  command,  electrical  power  system, 
and  sensor  subsystems.  The  paper  also  makes  the  distinction  between  satellites  as 
either  being  a  sensor-type  satellite  or  a  communications  satellite;  sensor-type  having 
a  sensor  based  payload  with  a  limited  communications  subsystem  and  communica¬ 
tion  having  no  sensors,  with  all  communications  payload.  Tactics  and  strategies  are 
related  to  parameters  “that  describe  the  size,  shape,  and  spatial  orientation  of  the 
satellite  orbit”  [17:501].  Strategies  consider  schedule  maintenance,  reactive  repair 
policies,  and  policies  involving  repairs  before  a  failure  is  anticipated.  Tactics  include 
’’expendable  and  reusable  transfer  vehicles”  and  ’’one  or  two  orbital  platforms  per 
scenario  to  serve  as  warehouses  and  transportation  nodes”  [17:501]. 

The  truncation  of  a  satellite’s  lifetime  is  stated  as  being  relatively  predictable, 
since  it  is  based  on  deterministic  factors  such  as  “mechanical  wearout  or  the  ex¬ 
haustion  of  consumables”  [17:502],  Also,  support  cost  was  determined  to  be  statis¬ 
tically  significant  with  the  following  four  constellation  parameters  and  five  satellite 
parameters:  number  of  satellites  in  constellation,  constellation  maintenance  time, 
cost  of  transportation  to  constellation,  location,  constellation  transportation  effi¬ 
ciency,  satellite  mass,  modularization  mass  penalty,  reliability,  truncation  lifetime, 
and  measure  of  value  of  retrieving  failed  satellites  and  modules.  They  used  fractional 
factorial  designs  to  minimize  run  lengths. 
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An  article  by  Hopp  [25]  discusses  incorporating  technological  improvements  in 
a  replenishment  model  for  a  general  system  by  basing  it  on  the  obsolescence  rather 
than  the  deterioration  of  the  equipment.  The  paper  states  that  past  replacement 
methods  are  based  on  equipment  wear  and  assume  that  replacements  are  identical 
to  the  previous  equipment.  Previous  papers  attempted  to  incorporate  the  obsoles¬ 
cence  of  systems  by  assuming  that  improving  technology  and  deteriorating  equipment 
were  deterministic,  while  [21]  attempted  to  model  the  technology  improvement  as  a 
stochastic  process  with  deterministic  equipment  deterioration.  Hopp  takes  the  next 
step  of  modeling  both  stochastic  technological  improvement  and  deterioration. 

A  single  technological  improvement  is  assumed  to  occur  in  an  unknown  time, 
resulting  in  a  recursive,  non-time  homogeneous  optimal  value  function  that  is  solved 
by  restricting  the  function  to  an  initial  state  and  applying  theory  of  forecast  hori¬ 
zons.  The  function  makes  the  decision  of  keeping  or  replacing  the  equipment.  New 
technology-based  replacements  are  assumed  to  be  better  than  the  legacy  equipment. 
The  paper  goes  on  to  show  the  effects  of  technological  change  on  the  replacement 
policy,  finding  that  available  technology  fosters  early  replacement  while  the  possi¬ 
bility  of  technology  being  available  in  the  future  fosters  incentive  to  keep  current 
technology  longer.  Methods  are  discussed  regarding  using  the  recursion  function 
with  time  in  reverse  to  calculate  the  expected  cost  of  a  sure  event  of  technological 
improvement  in  the  future,  given  the  use  of  old  technology  in  the  present. 

The  paper  continues  with  methods  of  computing  the  optimal  replacement  pol¬ 
icy  based  on  minimizing  the  net  present  cost  based  on  the  recursion  formulas.  They 
establish  that  a  forecast  horizon  will  exist  if  the  optimal  replacement  policy  is  unique, 
and  they  define  an  efficient  stopping  criteria  for  the  recursive  algorithm  that  includes 
the  chance  of  the  function  converging  on  a  solution  where  both  the  keep  and  replace 
portions  of  the  function  are  equal.  The  two  stopping  rules  Bound  Based  and  Bes 
and  Lasserre’s  rule  are  compared,  resulting  in  the  preference  to  the  Bound  Based 
rule  due  to  its  computational  efficiency  and  minimal  data  requirements. 
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2.f  Satellite  Capability  Measures 


Many  terms  have  been  used  to  capture  the  aspects  of  a  system  not  accounted 
for  by  its  reliability.  System  effectiveness  is  a  term  used  by  Ebcling  [14]  to  describe 
a  system’s  probability  of  performing  its  intended  mission,  incorporating  the  sys¬ 
tem’s  performance  and  availability.  Performability  is  a  related  measure  to  system 
effectiveness  which  is  a  composite  measure  of  the  performance  and  reliability  of  a 
system.  A  composite  measure  of  availability  and  performability  captures  a  system’s 
effectiveness. 

Meyer  [32]  gives  a  formal  presentation  of  the  ideas  of  performability,  effective¬ 
ness,  and  capability.  The  author  describes  perfomability  in  terms  of  a  probability 
measure  mapping  a  set  of  events  to  a  probability  space.  Performability  is  then  a 
function  mapping  an  accomplishment  set,  referring  to  measures  of  system  accom¬ 
plishment,  to  a  probability  measure.  He  defines  a  capability  function  as  the  mapping 
from  a  state  trajectory  from  a  stochastic  model  to  a  level  of  accomplishment. 

There  are  many  examples  of  papers  that  focus  on  measuring  qualities  of  satel¬ 
lite  constellations,  whether  they  be  termed  measures  of  effectiveness,  capability,  or 
performance.  In  [34],  Smith  models  the  performability  of  a  multiprocessor  system 
using  a  Markov  reward  model  (MRM),  which  incorporates  rewards  with  the  state- 
spaces  of  a  stochastic  process. 

The  thesis  by  Wilson  [42]  focuses  on  the  performability  of  a  geosynchronous 
weather  satellite  system,  although  he  refers  to  the  resultant  measure  as  a  measure 
of  effectiveness  (MOE).  He  uses  decision  analysis  methods  in  conjunction  with  pub¬ 
lished  sources  along  with  subject  matter  expert  opinions  on  what  defines  the  quality 
of  the  system  to  score  the  constellation.  A  list  of  attributes  are  first  defined  which 
represent  system  parameters  that  have  a  direct  impact  on  satellite  weather  systems. 
Then  these  parameters  are  evaluated  in  reference  to  area  coverage  of  the  earth.  Wil- 


2-8 


son  gives  a  hierarchical  model  of  the  weather  system’s  mission  which  ties  the  broad 
goals  of  the  mission  to  measurable  performance  parameters. 

A  thesis  by  Staats  [35]  developed  a  measure  of  a  satellite’s  utility  through 
multi-attribute  utility  theory.  The  reliability  of  the  satellite  is  incorporated  into  the 
measure  by  introducing  an  attribute  reflecting  the  SME’s  opinion  of  the  satellite’s 
expected  remaining  lifetime.  Although  the  reliability  is  not  explicitly  computed,  the 
resulting  measure  comes  closer  to  a  desired  measure  of  satellite  performance.  The 
terms  value  and  utility  are  often  used  interchangeably;  however,  this  thesis  addresses 
strictly  the  value  of  the  satellite  with  respect  to  its  contribution  to  the  constellation 
mission. 

2.5  Summary  and  Contributions 

Constellation-specific  methodologies  incorporating  maintenance  aspects  of  satel¬ 
lites  are  either  replenishment  policies  determining  the  level  of  spares  to  maintain,  or 
policies  deciding  when  to  optimally  replace  or  repair  a  satellite  based  on  net  present 
cost.  If  a  replenishment  policy  is  to  be  based  on  capability,  a  constellation  capability 
measure  will  be  required. 

Multi-attribute  utility  theory  has  been  the  most  common  method  of  computing 
measures  of  satellite  performance.  One  method  reviewed  attempted  to  incorporate 
reliability  into  the  satellite  value  using  utility  theory,  while  another  method  ignored 
aspects  of  reliability.  A  transient  model  of  satellite  degradation  was  not  found  in  the 
literature,  nor  was  a  measure  of  constellation  capability  based  on  the  availability  of 
satellite  functions. 

It  is  clear  from  the  literature  that  constellation  capability,  a  measure  of  the 
per  for  inability  of  the  satellite  constellation,  should  incorporate  both  aspects  of  its 
reliability  and  mission  effectiveness.  This  research  offers  a  stochastic  model  of  the 
degradation  process  of  a  constellation’s  satellites,  from  which  a  measure  of  availabil- 
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ity  is  computed.  Value  scores  of  satellite  capability  are  incorporated  into  the  measure 
using  value  theory  to  yield  a  reliability-based  performance  measure  of  constellation 
capability.  The  next  chapter  will  present  a  formal  mathematical  model  that  may  be 
used  to  compute  a  new  constellation  functional  capability  measure. 
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3.  Formal  Mathematical  Model 

This  chapter  will  give  a  description  of  the  formal  mathematical  model  for  con¬ 
stellation  capability.  By  modeling  the  status  of  each  satellite  in  the  constellation  as 
a  time-homogeneous,  continuous-time  Markov  chain  (CTMC),  and  by  using  decision 
analysis  techniques  to  subjectively  score  each  satellite,  we  derive  a  methodology  that 
provides  a  continuous  measure  of  capability  for  a  constellation.  Each  satellite’s  in¬ 
stantaneous  availability  and  value  score  are  evaluated  as  inputs  into  the  functional 
capability  measure.  The  instantaneous  availability  is  a  function  of  time  that  mea¬ 
sures  the  satellite’s  probability  of  being  available  at  a  future  point  in  time.  The  value 
score  is  a  measure  between  0  and  1  representing  the  satellite’s  capability  to  perform 
its  mission  within  the  constellation  based  on  a  subject  matter  expert’s  opinion. 

3.1  Constellation  Functional  Capability 

A  satellite  constellation  operates  in  the  environment  of  space  which  imposes 
stresses  on  its  satellites  that  may  potentially  result  in  a  degradation  of  functional 
capability.  It  is  impossible  to  completely  determine  how  long  a  satellite  constellation 
will  remain  in  a  functionally  capable  state  because  both  the  constellation  and  the  en¬ 
vironment  in  which  it  operates  in  are  not  completely  determined;  there  exist  random 
elements  in  the  system  and  in  the  environment  that  make  it  difficult  to  characterize 
the  behavior  of  the  system  over  time.  By  probabilistically  modeling  a  satellite  in  a 
constellation,  and  by  making  a  few  reasonable  assumptions  about  the  satellites  in 
the  constellation,  we  are  able  to  develop  a  mathematically  tractable  model  which 
gives  the  instantaneous  availability  and  a  measure  of  capability  for  the  satellite  con¬ 
stellation.  Instantaneous  availability  is  chosen  because  it  is  a  more  accurate  measure 
of  availability  for  systems  which  have  not  been  in  service  for  a  long  period  of  time. 

This  section  describes  the  constellation  functional  capability  measure 
composed  of  each  satellite’s  instantaneous  availability  score  and  value  score.  The 
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constellation  functional  capability  is  the  mathematical  average  of  the  functional  ca¬ 
pability  scores  for  the  satellites  within  the  constellation.  The  functional  capability 
measure  for  an  arbitrary  satellite,  namely  satellite  k,  is  the  product  of  its  instan¬ 
taneous  availability  at  time  t,  denoted  A ^  (t  —  t^^j,  and  its  value  score,  denoted 
V^k\t).  The  value  A®  (t  —  t^^j,  which  is  based  on  the  time  of  its  initialization  t^\ 
is  a  probability  and  as  such,  may  assume  a  value  between  0  and  1.  The  initialization 
time  tll'>  is  the  time  at  which  the  satellite’s  life  in  the  constellation  begins.  The 
value  score  V^k\t)  is  a  score  which  describes  the  functional  value  of  the  satellite  to 
the  constellation  mission  and  is  derived  through  multi-attribute  value  theory.  It  is 
dependent  on  the  time  the  constellation  is  evaluated  and  can  also  take  on  a  value 
between  0  and  1.  Given  the  satellite  constellation  consists  of  K  satellites, 

*(*)  =  v  X>w  (‘ -  4l))  r<l)«.  (3.i) 

Y  k= 1 

This  average  of  constellation  functional  capability  is  dependent  on  the  time  at 
which  the  “measurement”  is  taken,  and  on  the  set  of  criteria  comprising  the  value 
score  V^k\t).  The  capability  at  time  t,  <h(f),  will  be  a  measure  between  0  and 
1,  where  0  corresponds  to  a  constellation  that  has  no  functional  capability  and  1 
corresponds  to  a  constellation  that  has  full  functional  capability.  For  example,  if 
each  satellite  in  the  constellation  is  completely  available  at  time  t,  and  each  satellite 
is  also  scored  as  having  perfect  value,  then 

A {k)  (t  -  4fc))  =  1,  V{k\t)  =  1,  k  =  1,  2, . . . ,  K 

and 

m  =  (« -#>)  r<‘>(«)  =  ^E(i)(i)  =  jrK=1- 

k= 1  k= 1  ' 
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If  each  satellite  is  unavailable, 


A(fc)  (t  -  4fc))  =  0,  k  =  1,2,...,  K 

and 

m  =  ^£o-v«\t)  =  o. 

k= 1 

The  constellation  functional  capability  measure  <h(f)  captures  both  the  mission  con¬ 
tribution  and  the  availability  of  the  satellite,  where  the  score  takes  into  account 
aspects  of  each  satellites  reliability  and  fulfillment  of  operational  requirements. 

The  following  sections  will  focus  on  the  satellite  level  of  the  constellation  and 
will  discuss  the  assumptions  and  methodology  for  obtaining  the  instantaneous  avail¬ 
ability  and  value  score  for  the  satellites  when  the  degradation  status  of  each  satellite 
is  modeled  as  a  continuous-time  Markov  chain  (CTMC).  Since  the  literature  is  quite 
extensive  in  the  area  of  assigning  a  value  score  to  a  satellite,  the  methodology  for 
producing  a  value  score  is  covered  with  less  emphasis  than  the  methodology  for 
obtaining  its  instantaneous  availability  measure. 

3.2  Instantaneous  Availability 

A  model  for  instantaneous  availability  is  given  for  a  general  satellite  constel¬ 
lation  consisting  of  K  satellites,  where  the  degradation  state  of  a  satellite  can  be 
modeled  as  a  time-homogeneous,  continuous-time  Markov  chain  (CTMC);  the  failure 
and  repair  rate  distributions  of  the  satellite  states  are  restricted  to  exponential  dis¬ 
tributions.  Each  satellite  is  considered  independent  of  the  others  and  has  M ^  main 
functions  that  can  be  in  only  one  of  two  states:  operational  and  non-operational. 
A  more  general  model  is  presented  in  Chapter  4  for  which  the  function  failure  and 
repair  distributions  are  not  restricted  to  exponential  distributions. 
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3.2.1  Stochastic  Degradation  Model 

The  following  assumptions  apply  to  the  CTMC  model  for  computing  instan¬ 
taneous  satellite  availability: 

1.  The  constellation  under  consideration  contains  a  finite  number  (K)  of  inde¬ 
pendent  satellites,  each  with  a  set  of  finite  functions. 

2.  The  failure  and  repair  rates  of  each  function  of  the  satellite  are  assumed  to  be 
known. 

3.  The  distribution  of  failure  and  repair  times  is  exponential. 

4.  Multiple  satellite  function  failures  or  repairs  cannot  occur  in  any  instant  of 
time. 

5.  The  initial  state  of  the  satellite  constellation  is  assumed  to  be  known;  the 
states  the  satellite  may  be  in  are  assumed  to  be  definable  and  the  probabilities 
of  being  in  these  states  at  time  t^\  the  beginning  of  the  satellite’s  lifetime,  are 
assumed  to  be  known. 

6.  The  states  of  availability  of  the  satellite  can  be  determined  based  on  the  avail¬ 
ability  of  specific  functions. 

The  models  presented  in  this  chapter  and  in  Chapter  4  will  not  consider  the  pos¬ 
sibility  of  on-orbit  spares.  A  different  methodology  tailored  specifically  for  handling 
spares  and  for  addressing  long-term  constellation  replenishment  policies  is  provided, 
for  example  in  [8] ,  and  may  be  more  appropriate  for  assessing  constellation  availabil¬ 
ity  when  considering  satellite  replacement.  The  methodology  presented  in  this  and 
the  following  chapter  measures  the  constellation’s  current  capability  by  looking  at 
the  current  availability  levels  of  its  satellites.  Repairs  are  considered  at  the  satellite 
functional  level,  where  only  one  of  the  satellite’s  multiple  functions  may  fail  or  be 
repaired  at  any  instant  in  time.  In  order  to  incorporate  the  idea  of  spare  satellites 
into  the  methodology,  the  model  would  have  to  be  changed  to  allow  simultaneous 
failures  and  repairs  to  occur. 
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Figure  3.1  Sample  path  of  a  discrete  state  space  CTMC. 


In  order  to  obtain  an  instantaneous  availability  measure  for  a  satellite  constella¬ 
tion,  we  must  describe  the  constellation  in  terms  of  a  stochastic  model.  A  stochastic 
process  is  defined  by  Kulkarni  as  “a  probabilistic  model  of  a  system  that  evolves 
randomly”  [30:2],  A  stochastic  process  may  be  described  as  a  sequence  of  random 
variables  {X(t)  :  t  >  0}  taking  on  values  from  a  sample  space  S  with  X(t)  G  S 
for  each  t  >  0.  The  number  of  elements  in  set  S  is  referred  to  as  the  cardinality  of 
S  and  is  denoted  by  | S\.  The  cardinality  of  S  may  be  finite  or  infinite.  Stochastic 
processes  may  assume  discrete  or  continuous  values  over  discrete  or  continuous  time 
intervals.  For  certain  types  of  stochastic  processes,  information  such  as  the  proba¬ 
bility  the  system  will  be  in  a  particular  state  in  S'  at  a  specific  time  t,  or  the  long 
term  proportion  of  time  it  spends  in  a  given  state,  can  be  obtained  from  the  system. 

The  probability  distributions  of  the  duration  of  time  the  process  remains  in  its 
present  state,  referred  to  as  the  sojourn  time,  and  the  dynamics  of  state  transitions, 
distinguish  certain  types  of  stochastic  processes  from  others.  A  continuous-time 
Markov  chain  (CTMC),  denoted  {X(t)  :  t  >  0},  is  a  specific  stochastic  process 
which  describes  the  random  evolution  of  a  system  over  continuous  time  having  a 
countable  state  space. 

Figure  3.1  shows  a  graphical  representation  of  a  sample  path  of  a  CTMC, 
where  iq  is  the  sojourn  time  in  state  i.  In  the  figure,  the  system  starts  out  in  state  2 
(X  (0)  =  2)  and  the  process  remains  in  that  state  for  a  duration  of  T\  time  units  until 
the  process  makes  an  instantaneous  transition  to  state  3  at  time  S±  (A"(S'1h)  =  3). 
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The  CTMC  process  possesses  the  Markov  property  at  each  point  of  transition, 


where 


P{X(Sn+1)  =  j\X(Sn )  =  inXiSn-i)  =  in- 1, . . .  ,X(0)  =  i} 

=  P{X(Sn+1)=j\X(Sn)=in}. 

A  stochastic  process  is  Markovian  if  the  probability  of  transitioning  to  a  future  state 
only  depends  on  information  from  the  current  state.  The  process  remains  in  a  given 
state  for  a  random  amount  of  time,  which  is  distributed  exponentially.  The  CTMC 
is  time-homogeneous  if 


P{X(Sn+1)  =  j\ X(Sn)  =  *n,X(Sn_!)  =  in_!, . . .  ,X(0)  =  1} 
=  P{X(Sn+1)=j\X(0)=i}. 


This  implies  that  all  that  is  needed  to  find  the  probability  the  process  transitions  to 
state  j  is  to  know  the  initial  state  %  the  process  started  in.  Further  information  on 
CTMCs  is  given  in  Kulkarni  [30]. 

In  this  work,  we  are  modelling  the  stochastic  degradation  of  a  satellite  of  a 
constellation,  which  randomly  evolves  based  on  the  availability  of  its  functions. 
Define  Xm\t )  as  a  random  variable  which  describes  the  state  of  function  m  for 
satellite  k  at  time  t,  m  =  1,2, ,  M^k\  k  =  1,2, ... ,  K,  where  t  is  the  time  of 
inspection.  For  the  development  of  the  constellation’s  instantaneous  availability,  the 
superscript  index  will  indicate  which  satellite  in  the  constellation  we  are  referring  to 
and  the  subscript  index  will  indicate  which  function  of  the  satellite  we  are  referring 
to.  The  random  variable  Xm\t )  is  defined  as 


1,  if  function  m  is  available  on  satellite  k  at  time  t 
0,  if  function  m  is  not  available  on  satellite  k  at  time  t 


(3.2) 
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In  the  real  world,  an  unavailable  function  is  not  necessarily  a  failed  function  but 
for  this  methodology,  we  will  assume  that  a  function  that  is  in  a  failed  state  is 
synonymous  with  a  function  that  is  in  an  unavailable  state. 

For  modeling  the  state  of  the  satellite,  the  states  of  each  function  of  the  satellite 
can  be  grouped  into  a  row  vector  which  represents  the  state  of  satellite  k  at 

time  t: 

x(fc)(f)  =  x[  k\t)  xf\t)  ...  X^\t)  ■ 

Let  S  be  defined  as  the  state  space  of  a  satellite,  a  finite  set  containing  the 
possible  states  the  satellite  may  assume  at  any  point  in  time.  The  random  vector 
X®(£)  will  take  on  a  value  from  set  S  for  each  t  >  0.  The  cardinality  of  S  is 
|  S' |  =  2M{k\  since  this  represents  a  Bernoulli  trial  for  each  of  the  M ^  distinct 
satellite  functions. 

Each  satellite  is  assumed  to  be  observable  over  a  continuous  interval  of  time, 
where  it  may  only  transition  between  two  states  in  S  in  any  instant  of  time.  With 
this  model,  multiple  functions  may  not  fail  at  the  same  time;  for  repairs,  multiple 
functions  cannot  transition  from  a  failed  status  to  an  available  status  at  the  same 
time.  This  model  does  not  imply  that  repairs  cannot  be  made  concurrently,  only 
that  repairs  cannot  be  completed  simultaneously.  These  model  characteristics  allow 
the  satellite  to  be  modeled  as  a  CTMC  process,  denoted  {X^(£),£  >  0}. 

Transitions  from  one  state  in  S  to  another  are  based  on  the  discrete  event  of 
a  satellite  either  having  an  available  function  fail  or  having  an  unavailable  function 
repaired.  For  example,  a  satellite  with  M  =  3  functions  will  have  a  state  space  with 
23  =  8  elements  given  by 

(1,1,1), (1,0,1), (0,0,1), (1,0,0), 

(0,1,1),  (1,1,0),  (0,1,0),  (0,0,0) 

Table  3.1  shows  an  assignment  of  an  index  to  each  element  in  this  state  space. 
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Figure  3.2  Transition  rate  diagram  for  a  two-state  CTMC. 


Table  3.1  Index  assignment  for  a  three  function  satellite. 


State 

Combination 

Status 

1 

(1,1,1) 

All  functions  available 

2 

(o,i,i) 

Function  1  unavailable 

3 

(1,0,1) 

Function  2  unavailable 

4 

(1,1,0) 

Function  3  unavailable 

5 

(0,0,1) 

Functions  1  and  2  unavailable 

6 

(0,1,0) 

Functions  1  and  3  unavailable 

7 

(1,0,0) 

Functions  2  and  3  unavailable 

8 

(0,0,0) 

All  functions  unavailable 

A  CTMC  can  also  be  represented  graphically  with  a  transition  rate  diagram, 
a  directed  graph  where  each  node  of  the  graph  represents  an  element  of  S  and  each 
directed  arc  shows  a  possible  transition  that  may  occur.  Each  directed  arc  has  a 
rate  ql3  corresponding  to  the  rate  at  which  the  process  transitions  from  state  %  to 
state  j.  Figure  3.2  shows  an  example  of  a  transition  rate  diagram  for  a  two-state 
CTMC.  This  two-state  CTMC  is  often  referred  to  as  an  up-down  machine.  Figure 
3.3  shows  the  transition  rate  diagram  corresponding  to  the  index  assignment  from 
Table  3.1  along  with  all  of  the  possible  transitions  that  could  occur  between  the 
satellite  states. 

The  transition  probability  matrix,  denoted  P(f),  is  a  matrix  of  probabilities 
associated  with  the  CTMC  where  ptJ  (t)  is  the  {i ,  j)  entry  of  matrix  P (t) .  The  element 
Pij(t)  indicates  the  probability  of  the  process  transitioning  from  state  i  to  state  j  in 
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Figure  3.3  Transition  rate  diagram  for  a  three-function  satellite. 
t  time  units,  and  is  defined  as 

Pij(t)  =  P{x(t)  =  j\x(o)  =  i\,  i,jes 

for  a  time-homogenous  CTMC.  We  seek  the  probability  distribution  associated  with 
the  stochastic  process  {X(t),t  >  0}  being  in  each  of  the  states  of  the  sample  space. 
Let  aij(t)  denote  the  probability  that  the  process  is  in  state  j  at  time  t,  defined  as 

aj(t)  =  p{X(t)  =j},  j  e  S. 

The  probability  distribution  for  the  CTMC  ,  denoted  a(t),  is  a  row  vector  of  dimen¬ 
sion  |  (S' |,  where 

a(t)  =  ai(t)  a2(t )  ...  a\S\(t)  ■ 

The  row  vector  a(0)  is  the  initial  probability  distribution  which  corresponds  to  the 
probabilities  associated  with  the  process  being  in  each  of  the  possible  states  in  S  at 
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time  t  —  0,  where 


a(0)  —  «i(0)  a2(0)  ...  «|5|(0)  • 

An  important  theorem  in  [30]  and  [23]  states  that  a  time-homogeneous  CTMC 
is  completely  described  by  its  initial  probability  distribution  and  its  transition  prob¬ 
ability  matrix.  The  probability  mass  function  ct(t)  for  a  time-homogeneous  CTMC 
is  obtained  by  conditioning  on  the  probability  that  the  process  is  in  state  i  at  time 
t  =  0,  where  the  elements  are  computed  from 

aj(t)  =  ]T  P[X(t)  =  j\X(0)  =  i }  •  P{*(0)  =  i}, 

i£S 

resulting  in  the  probability  distribution 

a(t)  =  a(0)  •  P (t).  (3.3) 

Let  P^(f)  denote  the  transition  probability  matrix  for  satellite  k  where  p\j\t) 
is  the  (i,j)  entry  of  matrix  P (k\t)  and  indicates  the  probability  of  satellite  k  tran¬ 
sitioning  from  state  i  to  state  j  in  t  time  units.  The  elements  of  P(k\t)  are  defined 
as 

Pij\t)  =  P{x'k\t)  =  j|x<‘>(0)  =  o  i,j  €  s, 

since  the  satellite  is  modeled  as  a  time-homogenous  CTMC.  Let  cdfc)(f)  denote  the 
probability  distribution  of  the  status  of  satellite  k  at  time  t  and  let  cx^k\t^)  be 
the  initial  probability  distribution  which  contains  the  probabilities  associated  with 
satellite  k  being  in  each  of  the  possible  states  in  S  at  time  t  = 


where 

af]  =  P{X.(k\t(0k) )  =i}  ieS. 
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For  example,  for  satellite  k  =  4  with  M  ^  =  3  functions, 


0.5  0.5  0  0  0  0  0  0 


describes  satellite  4  starting  at  time  having  an  equal  probability  of  being  in  states 
1  and  2  and  having  no  probability  of  being  in  states  3  through  8.  If  a  checkout  of 
satellite  k  at  time  reveals  that  all  of  its  functions  are  available,  then  a^(i^)  will 
result  in  a  row  vector  with  the  first  element  being  equal  to  one  and  the  remaining 
elements  being  equal  to  zero. 

The  probability  mass  function  of  {X.(k\t),t  >  0}  is  then  given  by 

Eqn.  (3.3)  by  conditioning  on  the  probability  that  satellite  k  is  in  state  i  at  time 

,  ,(fc) 

t  =  t0  ,  given  as 

a(k\t)  =  c*(fc)(4fc))  •  P {k\t).  (3.4) 


Let  Q  denote  a  square  matrix  of  dimension  |  S\  which  contains  the  transition 
rates  from  one  state  to  another  within  set  S ;  qV)  is  the  (i,j)  entry  of  matrix  Q  and 
represents  the  rate  of  transition  from  state  i  to  state  j  within  the  (S)  states  (i  ^  j). 
The  matrix  Q  is  called  the  infinitesimal  generator  matrix.  The  elements  of  Q  are 
defined  by 

>=j  (3,5) 

Qij,  i  j 

where  the  ith  diagonal  entry  is  denoted  qt.  The  negative  of  the  diagonal  elements  of 
Q  are  the  rate  parameters  for  the  sojourn  time  distributions.  Since  CTMCs  in  this 
work  are  assumed  to  be  time-homogenous,  the  transition  rates  do  not  change  over 
time.  An  important  property  of  the  infinitesimal  generator  matrix  cited  in  [30]  is 


X^i  =  0-  (3-6) 

ieS 
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The  transition  rates  are  the  same  rates  listed  in  the  transition  rate  diagram  on 
each  arc.  In  Figure  3.2,  the  sojourn  time  associated  with  the  process  remaining  in 
state  1  is  exponentially  distributed  with  rate  A.  Similarly,  the  time  the  process  is  in 
state  2  is  exponentially  distributed  with  rate  //.  The  infinitesimal  generator  matrix 
corresponding  to  Figure  3.2  is 


Q 


-A  A 

H  —fi 


and  the  infinitesimal  generator  matrix  corresponding  to  the  transition  rate  diagram 
in  Figure  3.3  is 


—  (Ai  +A2  +  A3) 

Ai 

A2 

*3 

0 

0 

0 

0 

Ml 

-(Ml  +A2  +  A3) 

0 

0 

A2 

A3 

0 

0 

M2 

0 

-(M2  +  Ai  +  A3) 

0 

Al 

0 

A3 

0 

M3 

0 

0 

-(M3  +  Ai  +^2) 

0 

Ai 

A2 

0 

0 

M2 

Ml 

0 

-(mi  +M2  +  A3) 

0 

0 

A3 

0 

M3 

0 

Ml 

0 

-(mi  +M3  +A2) 

0 

A2 

0 

0 

M3 

M2 

0 

0 

-(M2+M3  +  Ai) 

Ai 

0 

0 

0 

0 

M3 

M2 

Ml 

-(mi  +M2  +M3) 

where  A*  is  the  failure  rate  of  function  i,  and  yU*  is  the  repair  rate  of  function  i. 
Recall  that  the  negative  of  the  diagonal  elements  of  Q  are  the  rate  parameters  for 
the  sojourn  time  distributions. 

The  transition  rate  matrix  and  the  probability  transition  matrix  are  related 
by  the  following  equations,  also  known  as  the  forward  and  backward  Chapman- 
Kolmogorov  equations 

jP  (t)  =  P(t)Q  =  QP(t).  (3.7) 

Eqn.  (3.7)  is  a  system  of  homogeneous,  first  order  linear  differential  equations  with 
initial  condition 

P(0)  =  I,  (3.8) 

where  I  is  the  identity  matrix  with  dimension  | S\.  By  definition  of  a  CTMC,  at 
time  t  =  0  it  is  a  sure  event  that  the  system  will  either  be  in  state  1  or  state  2;  the 
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process  cannot  be  in  two  states  at  the  same  time.  An  exact  solution  to  Eqn.  (3.7) 
can  be  obtained  using  methods  for  solving  first  order  differential  equations.  Once  the 
transition  probability  matrix  is  obtained,  we  can  solve  for  the  probability  distribution 
associated  with  the  process. 

Let  S'  be  the  subset  of  S  which  describes  the  states  of  the  satellite  when  it  is 
considered  to  be  available  to  perform  its  mission,  defined  as 


S'  —  { <Ti ,  CT2 ,  .  .  .  ,  <T|S'|}  C  S, 


where  <7*  refers  to  the  ith  element  of  S'  according  to  its  index  assignment.  The  in¬ 
stantaneous  availability  A ^  (t  —  of  satellite  k  initialized  at  time  t^\  describes 
the  probability  of  the  satellite  being  in  a  state  of  availability  at  time  f,  which  is 
equivalent  to  the  sum  of  the  probabilities  of  the  state  of  satellite  k  being  equal  to 
the  elements  in  set  S' ; 

1  (f  -  =  V  a,(t  -  «<*>)  =  Y,  P(X{t\t  -  i<*>)  =  <r,).  (3.9) 

jeS'  i=i 

3.2.2  Solution  Methods 

A  review  of  two  methods  for  solving  a  system  of  linear  differential  equations 
for  the  probability  transition  matrix  is  given  below,  one  using  direct  integration, 
and  the  other  using  Laplace  transforms.  An  example  of  solving  a  two-state  system 
follows,  and  [16]  provides  further  explanation. 

The  first  method  discusses  using  direct  integration  to  solve  for  the  probability 
transition  matrix  given  the  rate  matrix.  We  start  with  the  Kolmogorov  forward 
equations  in  standard  form  where 

^P(f)  +  -QP(f)=0.  (3.10) 
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The  integrating  factor,  denoted  p(t),  can  be  obtained  by  solving  the  equation 


pit)  =  exp 


which  yields 


p{t)  =  exp  (— Q t) 


since  Q  is  a  constant.  Multiplying  both  sides  of  Eqn.  (3.10)  by  this  integrating  factor 
gives 


exp  (— Q t) 


jP (t)  -  QP (() 


0. 


(3.11) 


We  can  recognize  the  left  side  of  Eqn.  (3.11)  as  the  derivative  of  a  product  of  the 
integrating  factor  exp  (— Q t)  and  P(f).  Thus  Eqn.  (3.11)  becomes 


jf  {exp  (— Qt)P(t)}  =  0. 


When  we  integrate  both  sides  with  respect  to  t, 


exp  (— Qt)P(t)  =  C,  (3-12) 

where  C  is  a  constant,  square  matrix  of  dimension  \S\.  From  our  initial  condition, 
the  matrix  P(0)  =  I.  Using  this  fact,  we  can  solve  for  C: 

C  =  exp  (— Q  ■  0)P(0)  =  I. 

After  substituting  I  for  C  into  Eqn.  (3.12)  and  solving  for  P(t),  we  obtain 

P(t)  =  exp(Qt), 
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where  exp  (Q t)  is  the  matrix  exponentiation  of  Q t  and  is  defined  as 


exp(Qt)  =  f 

'  n\ 

n= 0 


i+E 

n=  1 


m 

n\ 


(3.13) 


The  use  of  Laplace  transforms  is  a  second  method  to  solve  Eqn.  (3.7)  for  P (t). 
This  involves  transforming  the  system  from  the  time  domain  to  the  complex  number 
domain,  where  the  system  of  ordinary  differential  equations  becomes  a  system  of 
linear  equations  that  is  much  easier  to  solve.  After  the  transform  solution  is  obtained 
in  the  transform  domain,  the  system  is  transformed  back  to  the  time  domain  through 
an  inverse  Laplace  transform.  Let  fit)  be  an  absolutely  integrable  function  on  the 
positive  real  line  and  let  f*(s)  =  £{f(t)}  denote  the  Laplace  transform  of  f(t)  given 
by 

/»oo 

f*(s )  =  /  e~stf{t)dt,  t  >  0,  s  G  C,  5R(s)  >  0,  (3.14) 

Jo 

where  C  denotes  the  set  of  complex  numbers  and  5R(s)  denotes  the  real  part  of  the 
complex  number  s.  The  Laplace  transform  of  the  probability  transition  matrix  P  (t) 
is  denoted  by  P*(s),  where  its  (i,j) th  element  p*-(s)  is  defined  as 


e  stpij(t)dt. 


We  will  apply  the  Laplace  transform  to  both  sides  of  Eqn.  (3.7).  From  prop¬ 
erties  of  Laplace  transforms  [30] , 


£{spw}  =  sp-w-p<°)- 


The  Laplace  transform  of  a  function  multiplied  by  a  constant  is  equal  to  the  constant 
multiplied  by  the  Laplace  transform  of  the  function.  In  the  case  of  Eqn.  (3.7),  Q  is 
a  constant: 

£{QP(f)}  =  Q£{P(t)}  =  QP*(s). 
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Applying  the  Laplace  transform  to  Eqn.  (3.7)  yields 


sP*(s)  -P(0)  =  QP*(s). 

Using  the  initial  condition  P(0)  =  I  with  factoring  yields 

(si  —  Q)P*(s)  =  I. 

After  left-multiplying  both  sides  of  the  equation  by  (si  —  Q)_1,  we  have 

P*(s)  =  (si  -  Q)"1.  (3.15) 

Eqn.  (3.15)  is  a  matrix  of  rational  functions  in  s.  The  inverse  Laplace  transform 
can  be  applied  to  P*(s)  to  obtain  an  exact  solution  via  properties  of  inverse  Laplace 
transforms  and  partial  fraction  decomposition  [30].  The  inverse  Laplace  transform 
of  P*(s),  denoted  £_1{P*(s)},  brings  the  probability  transition  matrix  back  into  the 
time  domain  where 


P(t)  =  £_1{P*(s)}  =  ^{(sl  -  Q)"1}.  (3.16) 

The  formal  equation  for  the  inverse  Laplace  transform,  also  known  as  the  Bromwich 
or  Fourier-Mellin  integral,  is  given  by 

-1  nc+ioo 

£_1{ /*(* )}  =  f(t)  =  —  /  estf*(s)ds,  t  >  0,  (3.17) 

Z7TZ  J c—ioo 

where  i  =  \f— T.  The  exact  solution  to  P(t)  can  be  obtained  by  observation,  by 
applying  some  basic  properties  of  inverse  Laplace  transforms  after  partial  fraction 
decomposition  has  been  performed. 

The  solution  obtained  using  an  integrating  factor  and  using  the  Laplace  trans¬ 
form  technique  yield  the  same  solution.  An  approximation  to  P(f)  can  be  obtained  by 
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applying  numerical  methods  to  solve  Eqn.  (3.13)  or  Eqn.  (3.15).  Many  matrix-based 
computer  software  packages  come  equipped  with  matrix  exponentiation  functions 
that  can  be  applied  to  approximate  the  solution  to  Eqn.  (3.13).  For  example,  Mat- 
Lab  release  13  comes  equipped  with  several  variations  of  a  function  that  performs 
matrix  exponentiation.  Numerical  algorithms  for  performing  the  inverse  Laplace 
transform  are  numerous  and  stem  from  a  seminal  paper  from  Dubner  and  Abate 

[13]- 

To  calculate  the  instantaneous  availability  of  a  system,  a  subset  must  be  de¬ 
fined  in  the  system’s  state  space  that  is  equivalent  to  the  system  being  available. 
After  this  subset  is  defined,  transient  analysis  can  be  performed  to  determine  the 
probabilities  associated  with  the  system  being  in  those  states  of  availability  at  time  t; 
the  instantaneous  availability,  A(t),  can  be  determined  by  summing  the  probabilities 
associated  with  those  states.  Let  S'  be  the  subset  of  S  which  describes  the  states  of 
availability: 

S'  =  {ai.a-2,..  ■ ,  cr|s,,| }  C  S. 

The  value  A(t)  is  then  equal  to  the  sum  of  the  probabilities  of  the  system  being  in 
set  S'  at  time  t: 

^(*)  =  =  I>(x(i)  =  <7,).  (3.18) 

jeS’  i= i 

Next,  we  discuss  the  CTMC  analysis  for  determining  the  probability  mass 
function  and  the  instantaneous  availability  for  the  two-state  CTMC  in  Figure  3.2, 
where  S  =  {1,2}.  For  calculating  the  probability  mass  function  a(t),  we  need  the 
initial  probability  distribution  vector,  the  infinitesimal  generator  matrix,  and  the 
desired  time  t.  For  calculating  the  instantaneous  availability  A(t),  we  need  to  know 
what  states  in  S  are  indicative  of  the  system  being  available. 
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X(t) 
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Figure  3.4  Sample  path  of  a  two-state  CTMC. 

Let  the  system  be  available  when  it  is  in  state  1  and  let  the  system  start  in 
state  1  at  time  t  =  0  by  having  an  initial  probability  distribution  vector 


q(0)  = 


1  0 


The  system  is  up  when  it  is  in  state  1  and  down  when  it  is  in  state  2.  Figure  3.4 
shows  a  sample  path  of  an  up  down  machine  starting  in  state  1.  Let  the  system  have 
exponential  distributions  for  its  time  in  each  state,  where  A  and  /i  refer  to  the  rates 
of  failure  and  repair,  respectively,  resulting  in  an  infinitesimal  generator  matrix 


,  d  ~d  J 

We  will  now  solve  Eqn.  (3.7)  for  P (t).  The  forward  Chapman  Kolmogorov  equation 

4p(i)  =  pWQ 

becomes 


d 

Pu(t)  Pl2(t) 

Pu(t)  Pn(t) 

-A  A 

dt 

P2l(t)  P22{t) 

P2l(t)  p22(t ) 

ji  —p 

— Apn(f)  +  ppi2(t)  Apn(f)  -  ppl2(t) 

-Ap2l(t)  +  PP22{t)  Xp2l(t)  -  pp22(t) 
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which  is  equivalent  to  the  linear  differential  equations 


Jfp  u{t)  =  -Xpn(t)  +  npi2(t) 

(3.20) 

^Pi2(t)  =  Xpu(t)  -  ppr2(t) 

(3.21) 

^P2l(t)  =  ~Xp2i(t)  +PP22(t) 

(3.22) 

JfP22(t)  =  Xp2l(t)  -  fip2z(t) 

(3.23) 

with  initial  conditions 

Pii(O)  =  p22(0)  =  1  and  pi2(0)  =  p2i(0)  =  0. 

From  an  axiom  of  probability,  the  rows  of  P(f)  have  to  sum  to  1: 

P11  +  P12  =  P21  +  P22  =  1-  (3.24) 

Applying  this  rule  of  complements  allows  us  to  solve  for  all  of  the  probability  values 
by  only  solving  one  probability  value  from  each  row  of  P(f).  We  will  show  the 
solution  for  pn  and  pi2. 

Using  substitution  from  Eqn.  (3.24)  and  arranging  in  standard  form,  Eqn.  (3.20) 
becomes 

^_Pn{t)  +  (X  + p)pn{t)  =  p.  (3.25) 

Our  integrating  factor  is  observed  to  be 

p(t)  =  exp  ((A  +  p)t), 

and  multiplying  both  sides  of  Eqn.  (3.25)  by  this  integrating  factor  yields 

exp  ((A  +  p)t)  ^pn(i)  +  (A  +  p)pu(t)  =  exp  ((A  +  p)t)p.  (3.26) 
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We  can  recognize  the  left  side  of  Eqn.  (3.26)  as  the  derivative  of  a  product  of  the 
integrating  factor  exp  ((A  +  p)t)  and  puit).  Thus  Eqn.  (3.26)  becomes 


d 

dt 


{exp  ((A  +  /t)f)pn(f)} 


exp  ((A  +  p)t)p 


and  integrating  both  sides  with  respect  to  t  gives 

exp  ((A  +  p)t)pu(t)  =  — ^ —  exp  ((A  +  p)t)  +  c.  (3.27) 

A  +  p 

From  our  initial  condition,  pn(0)  =  1  can  be  used  to  solve  Eqn.  (3.27)  for  c  resulting 
in 

A 

C  ~  T - • 

A  +  p 

After  substituting  the  result  for  c  into  Eqn.  (3.27)  and  solving  for  pn(t),  we  obtain 

Puit)  =  t-j—  +  exp  (-(A  +  p)t). 

A  +  p  A  +  p 

From  Eqn.  (3.24), 


Pi2(t)  =  l~Pn(t) 
p 


=  1  - 
A 


+ 


A 


■  exp  (-(A  +  p)t) 


A  T  /i  A  T  /i 

exp  (-(A  +  p)t). 


A  T  p  At  p 


In  a  similar  manner,  solutions  for  equations  (3.22) 
resulting  in  the  matrix 


and  (3.23)  can  be  obtained 


p« 


-JL.  + 


A  -\-fi 


a+7  exp 

W7exP 


(—(A  T  p)t ) 
(—(A  T  p)t) 


A 

A  -\-fi 


A+/i. 


T 


a+7  exp 

W7exP 


(-(A 

(-(A 


+  d)t) 
+ 
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The  probability  mass  function  can  be  obtained  by  multiplying  the  initial  prob¬ 
ability  distribution  to  P  (t),  resulting  in 


a(t)  =  a(0)P(t) 


1  0 


xfc  +  MT  exP  (-(A  +  /*)*) 
xfe  -  A^exp(-(A  +  ^)t) 


A 

\+u 


A 


The  availability  of  the  system  is  calculated  by  summing  all  of  the  probabilities 
of  the  states  which  indicate  the  system  is  available.  In  this  case,  the  system  is 
available  when  it  is  in  state  1.  Therefore, 

A(t)  =  ~T~7  +  VT-  exp  +  vW 

A  T  /i  A  T  /i 

is  the  probability  the  system  will  be  available  at  time  t  with  failure  and  repair  rates 
A  and  y,  respectively. 


3.3  Limiting  Availability 

Methods  for  obtaining  the  limiting  availability  for  a  system  modeled  as  a 
CTMC  are  given  below  for  comparison  purposes.  Chapter  5  will  present  numeric 
examples  with  calculations  of  both  the  limiting  availability  and  the  instantaneous 
availability  to  show  that  these  methods  produce  significantly  different  results  when 
assessing  a  constellation’s  functional  capability. 

As  was  previously  stated,  limiting  availability  is  a  long-run  average  of  the  pro¬ 
portion  of  time  a  system  is  available  over  its  entire  life  span.  Given  the  instantaneous 
availability  A(t),  the  limiting  availability  A  is  defined  as 

A  =  lim  A(t). 

t— XX) 
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In  terms  of  satellite  availability,  let  and  A^k\t)  denote  the  limiting  availability 
and  the  instantaneous  availability  for  satellite  k,  respectively.  Then 

\s'\ 

A{k)  =  lim  A{k\t)  =  lirn  V  P(X(fc)(f)  =  a). 

t — XX)  t—> oo  ‘  ^ 

i=  1 

The  methods  for  obtaining  the  limiting  availability  for  a  system  modeled  as  a 
CTMC  consist  of  obtaining  the  probability  distribution  for  the  satellite  when  it  has 
reached  a  limiting  state  of  equilibrium,  then  summing  up  the  long-term  probabilities 
associated  with  states  of  availability.  Let  p'^  denote  a  row  vector  consisting  of 
probabilities  pf}  defined  as 

p[k)=  ]im  P(XW(t)  =  ai),  a^eS,  i  -1.2, . , 

t— X) o 

(k) 

The  probability  p\  represents  the  long  term  probability  the  satellite  is  found  in 
state  i  and  can  also  be  view  as  the  proportion  of  time  the  system  spends  in  state  % 
over  its  entire  life.  The  system  of  equations 

p(fc)Q(fc)  =  0 

p(fc)e  =  1 

is  solved  for  p^,  where  is  the  rate  matrix  associated  with  satellite  k,  0  denotes 
a  row  vector  of  zeros  and  e  denotes  a  row  vector  of  ones.  The  first  equation  im¬ 
plies  the  condition  of  the  satellite  being  in  a  state  of  equilibrium  where  the  rate  of 
probability  flow  from  one  state  to  another  is  equal  to  zero.  The  second  equation  is 
the  normalizing  condition,  corresponding  to  the  axioms  of  probability.  The  limiting 
availability  can  then  be  obtained  by  summing  the  probabilities  associated  with  states 
of  availability  where 

\s’\ 

i= 1 
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3-4  Satellite  Value  Score 


This  section  gives  an  overview  of  one  possible  methodology  for  creating  the 
satellite  value  score  V^(t),  the  value  of  satellite  k  to  support  the  overall  constellation 
mission  at  time  t.  This  measure  is  included  to  account  for  scenarios  in  which  a 
satellite  may  be  functioning  at  its  specified  design  but  may  not  be  contributing  to 
the  current  mission  of  the  satellite  constellation.  This  situation  may  occur  dne  to 
problems  associated  with  the  acquisition  of  the  satellite,  such  as  a  lengthy  acquisition 
schedule  resulting  in  the  fielding  of  a  satellite  with  obsolete  technology  or  a  satellite 
that  no  longer  meets  mission  requirements  due  to  a  dynamic  threat  environment. 

For  this  particular  method,  the  satellite  value  score  is  obtained  through  the 
use  of  multi-attribute  value  theory,  which  describes  the  methodology  of  assigning 
numerical  value  to  a  subject  through  the  direct  evaluation  of  quantitative  attribute 
scores  or  through  the  scoring  of  qualitative  attributes  by  a  subject  matter  expert 
(SME).  In  this  application,  the  satellite  value  score  is  a  normalized  quantity  that 
represents  the  opinion  of  a  SME  on  the  operational  value  of  a  satellite,  where  the 
attributes  reflect  a  satellite’s  contribution  to  the  overall  constellation  mission. 

The  value  score  can  take  on  many  different  forms.  An  additive  form  for  this 
value  is  given  in  [7]  and  [29]  and  is  computed  for  each  satellite  as 

N  N 

V{k)(t)  =  ^2wf(t)vf)(yi),  ^w\k\t)  =  1,  (3.28) 

i=l  i—  1 

where  N  denotes  the  number  of  attributes  used  to  score  the  satellite,  w\k\t)  is  the 
normalized  weighting  factor  for  attribute  i  of  satellite  k  at  time  t,  and  v('k\yi)  is  the 
value  function  of  attribute  i  for  satellite  k  evaluated  with  score  ?/*.  In  order  to  use 
Eqn.  (3.28)  to  obtain  a  value  score,  the  attributes  are  assumed  to  meet  the  condi¬ 
tion  of  additive  independence  given  in  [7].  Obtaining  the  value  score  requires  first 
determining  which  attributes  to  include,  creating  scales  for  qualitative  attributes, 
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determining  an  appropriate  value  function  for  each  attribute  which  maps  the  at¬ 
tribute  scale  to  a  normalized  attribute  value,  then  rank-ordering  and  assigning  a 
weight  of  importance  to  each  attribute.  The  value  score  will  be  obtained  by  sum¬ 
ming  the  product  of  the  attribute  scores  and  the  attribute  weights.  The  input  of  a 
SME  familiar  with  both  the  operation  of  the  satellites  and  the  mission  that  depends 
on  the  constellation  will  be  critical  to  establishing  meaningful  satellite  value  scores. 

3-4-1  Attribute  Value  Functions 

The  first  step  of  this  methodology  involves  determining  a  set  of  applicable  at¬ 
tributes  that  reflects  the  satellite’s  ability  to  fulfill  current  mission  requirements.  An 
attribute  can  be  quantitative  or  qualitative.  For  attributes  of  a  quantitative  nature 
such  as  satellite  performance  metrics,  an  appropriate  range  must  be  determined  to 
create  the  attribute  scale.  For  attributes  of  a  qualitative  nature,  a  categorical  scale 
should  be  created  to  differentiate  between  different  levels  of  the  attribute,  where  a 
numeric  scale  can  be  created  from  the  worst  and  best  subjective  categories.  The 
attribute  scales  are  generated  to  create  a  domain  for  the  attribute  value  functions. 

For  example,  in  [35],  an  attribute  describing  a  satellite’s  contribution  to  the 
parent  constellation  mission  is  assigned  a  scale  from  0  to  1,  where  the  satellite  is 
given  an  attribute  score  of  0.0  if  it  is  a  spare  or  does  not  offer  any  support  to  the 
constellation  mission,  0.25  if  it  only  has  secondary  impact  on  the  constellation  mis¬ 
sion,  0.5  if  it  generally  impacts  the  constellation  mission,  0.75  if  it  strongly  impacts 
the  constellation  mission,  and  1.0  if  it  is  critical  to  the  constellation  mission.  In 
this  particular  case,  the  attribute  scale  is  normalized  but  this  is  not  a  necessary 
requirement. 

Next,  the  value  function,  denoted  v<jfk\yi),  is  obtained  for  attribute  i  of  satellite 
k.  The  value  function  is  created  with  the  help  of  a  SME  and  yields  a  normalized 
value  representing  the  satellite’s  value  relative  to  the  attribute.  After  value  functions 
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are  created  for  each  attribute,  the  satellite  attributes  can  be  scored  by  a  SME  using 
interview  questions,  where  yt  denotes  the  score  for  attribute  i. 

For  example,  given  attributes  with  normalized  scales,  a  value  function  of  an 
exponential  form  can  be  created  by  the  SME  determining  where  on  the  attribute 
scale  a  value  score  of  0.5  would  correspond  to.  Based  on  this  measure,  denoted  here 
as  y[,  vfc\y'i)  =  0.5.  The  exponential  value  function  would  then  be  given  as 


(Vi) 


1  -  exp(cf  }2/j) 
l-exp(cSfc))  ’ 


(3.29) 


(k) 

where  c\  ;  is  a  constant  that  can  be  computed  by  applying  numerical  techniques  such 
as  the  bisection  method  or  Newton’s  method  to  solve  the  equation 


1  -  exp(cf  }y') 
1  -  exp(cf }) 


(3.30) 


The  constant  is  referred  to  in  [35]  as  the  risk  attitude  constant  (RAC).  Regarding 
the  SME’s  chosen  value  for  y[,  the  value  function  is  linear  if  y[  =  0.5,  concave  if 
y\  <  0.5,  and  convex  if  y'i  >  o.5. 

For  attribute  scales  which  may  not  be  normalized,  a  value  function  of  the  form 


(Vi) 


1  H-  a 


-biVi 


(3.31) 


may  be  used,  where  a*  and  5*  are  constants  associated  with  attribute  i.  Eqn.  (3.31) 
is  referred  by  [11]  as  an  S-shaped  curve.  In  order  to  solve  for  these  constants,  more 
than  one  score  on  the  attribute  scale  representing  value  must  be  solicited  from  an 
SME  (i.e.,  multiple  scores  of  y[).  The  constants  a*  and  bi  can  then  be  solved  for 
using  numerical  methods. 
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3-4-2  Attribute  Weights 

Once  every  attribute  has  a  corresponding  value  function,  attribute  weights 
can  be  assigned,  allowing  preference  between  attributes  to  be  accounted  for  in  the 
satellite  value  score.  As  stated  earlier,  w\k\t)  refers  to  the  normalized  weighting 
factor  of  attribute  i  for  satellite  k  at  the  time  of  evaluation  of  the  constellation. 
This  measure  is  time-dependent  to  account  for  changes  in  the  order  of  preference  of 
attributes  for  constellation  missions  of  a  dynamic  nature. 

Normalized  attribute  weights  are  necessary  components  of  the  satellite  value 
score  and  methods  for  determining  these  weights  have  been  studied  extensively  in 
the  literature.  For  completeness,  we  point  out  three  methods  described  in  [7]  for 
assessing  attribute  weights:  the  pricing-out  method,  the  swing-weight  approach,  and 
the  lottery-comparison  technique.  A  more  mathematical  explanation  of  the  swing- 
weight  approach  and  the  pricing-out  method  is  given  in  [29],  where  attribute  weights 
are  referred  to  as  scaling  constants. 

3. 5  Model  Summary 

This  chapter  showed  the  formulation  of  the  constellation  functional  capability 
measure  by  first  introducing  the  equation  for  the  measure,  giving  a  review  of  CTMC 
analysis,  then  explaining  the  two  components  of  the  capability  measure  in  terms  of 
an  individual  satellite’s  instantaneous  availability  and  its  functional  value  score.  A 
model  was  given  for  satellite  instantaneous  availability  where  the  satellites  are  mod¬ 
eled  as  time-homogeneous  CTMCs  corresponding  to  exponential  failure  and  repair 
distributions.  Next,  a  method  for  computing  the  limiting  availability  was  given  for 
comparison  purposes.  Finally,  the  satellite  value  score  was  derived,  giving  an  ex¬ 
ample  of  one  methodology  that  may  be  used  to  assign  a  value  to  a  satellite,  based 
on  criteria  determined  by  a  subject  matter  expert,  which  reflects  the  capability  of 
a  satellite  fulfilling  a  specific  mission  requirement.  The  next  chapter  will  discuss 
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an  extension  to  the  methodology  presented  in  this  chapter  for  computing  the  in¬ 
stantaneous  capability  of  a  satellite  constellation  given  its  function  lifetimes  are  not 
exponentially  distributed. 
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4.  Semi-Markov  Degradation  Process 

This  chapter  describes  an  alternate  (and  more  general)  method  for  computing 
the  instantaneous  capability  of  a  constellation  containing  satellites  with  generally 
distributed  failure  and  repair  times.  The  methodology  is  similar  to  that  of  Chap¬ 
ter  3,  however,  the  satellite’s  instantaneous  availability  is  acquired  by  modeling  its 
degradation  status  as  a  time-homogenous  semi-Markov  process  (SMP)  where  the 
failure  and  repair  time  distributions  are  not  restricted  to  exponential  distributions. 
The  satellite  value  scores  are  computed  in  the  same  manner  as  described  in  the 
previous  chapter. 

4-1  Semi-Markov  Process  (SMP)  Model 

Recall  that  X^(t)  is  the  random  vector  describing  the  state  of  satellite  k  at 
time  t,  where  each  element  in  the  vector  indicates  the  state  of  availability  of  the 
specific  functions  of  the  satellite.  State  space  S  contains  all  of  the  combinations 
of  states  satellite  k  may  assume.  The  stochastic  process  (X^(f)  :  t  >  0}  evolves 
randomly  assuming  values  in  state  space  S. 

In  order  to  obtain  the  constellation’s  functional  capability  measure,  each  satel¬ 
lite’s  instantaneous  availability  must  be  obtained.  The  instantaneous  availability 
measure  for  a  satellite  constellation  that  has  non-exponential  failure  and  repair  dis¬ 
tributions  must  be  described  in  terms  of  a  semi-Markov  process  (SMP).  Modeling  the 
satellite  as  a  SMP  requires  knowledge  of  the  distributions  of  the  state  holding  times 
of  the  satellite  and  the  transition  probabilities  associated  with  each  of  the  states  in 
S.  The  following  assumptions  apply  to  the  SMP  model  of  satellite  degradation: 

1.  The  constellation  contains  a  finite  number  of  independent  satellites,  each  with 

a  finite  set  of  functions. 


4-1 


2.  The  transition  rates  between  states  of  availability  are  assumed  to  be  estimable 
from  observed  data. 

3.  The  probability  distributions  of  sojourn  durations  in  individual  states  are 
known  or  may  be  estimated  (parametrically). 

4.  Multiple  function  failures  or  repairs  cannot  occur  in  any  instant  of  time. 

5.  The  initial  state  of  the  satellite  constellation  is  assumed  to  be  known;  the 
states  associated  with  the  different  combinations  of  states  the  satellite  may 
be  in  are  assumed  to  be  definable  and  the  probabilities  associated  with  the 
satellite  being  in  these  states  when  the  satellite  first  came  on  line  are  assumed 
to  be  known. 

6.  The  availability  status  of  the  satellite  can  be  determined  based  on  the  avail¬ 
ability  of  specific  functions. 

A  short  review  of  SMPs  is  given  before  the  SMP  degradation  model  is  pre¬ 
sented.  A  description  of  renewal  processes  including  the  more  specific  Markov  re¬ 
newal  processes  is  first  given,  as  SMPs  are  defined  by  an  embedded  Markov  renewal 
process.  Further  information  on  SMPs  is  given  in  Kulkarni  [30]. 

A  renewal  process  is  a  counting  process  which  counts  the  number  of  events 
from  a  renewal  sequence.  Let  Sn  be  the  time  of  the  nth  event  and  let  N(t)  be 
the  number  of  events  up  to  time  t.  If  the  times  between  events  are  independent  and 
identically  distributed  random  variables,  then  {Sn  :  n  >  0}  is  a  renewal  sequence  and 
{N(t)  :  t  >  0}  is  a  renewal  process  [30].  Similarly,  a  Markov  renewal  process  counts 
the  number  of  events  from  a  Markov  renewal  sequence.  A  Markov  renewal  sequence 
is  a  sequence  of  bivariate  random  variables  {( Yn ,  Sn )  :  n  >  0}  where  Yn  indicates  the 
observation  of  the  process  at  the  nth  transition  epoch  Sn.  Let  { X(t )  :  t  >  0}  denote 
a  continuous-time  stochastic  process  where  X  ( t )  G  S,  and  let  Xn  denote  the  state  of 
the  process  after  the  nth  transition  epoch  Sn,  where  Xn  =  X(S^).  From  [30],  the 
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sequence  {(Xn,  Sn)  :  n  >  0}  is  a  Markov  renewal  sequence  if 

P{Xn+ 1  j .  >S,n_|_i  Sn  ^  t|An  f,  An_i, . . . ,  Ao,  Sn— i,  Sn—2i  ■  •  •  1 S0} 

=  P{Xn+1  =  j,  Sn+ 1  -  Sn  <  t\Xn  =  i}  (Markov  Property).  (4.1) 

The  Markov  renewal  sequence  is  time-ho7nogeneous  if 

P{Xn+1  =  J ,  Sn+i  -  Sn  <  t\Xn  =  i}  =  P{X1  =  J,  S,  <  t\X0  =  i}.  (4.2) 

The  Markov  property  implies  that  the  probability  that  the  process  transitions 
from  state  i  to  state  j  at  any  transition  epoch  can  be  determined  by  knowing  the 
current  state,  regardless  of  the  history,  along  with  the  sojourn  time  distribution  in 
state  i.  Time  homogeneity  implies  that  these  transition  probabilities  and  associated 
sojourn  time  distributions  do  not  change  with  time. 

A  continuous-time  stochastic  process  (A"(t)  :  t  >  0}  is  a  SMP  if  it  has  an 
embedded  Markov  renewal  sequence  {(An,^)  :  n  >  0}  where  Xn  is  defined  to  be 
equal  to  X(S'+),  the  state  after  the  nth  transition  epoch  Sn.  An  SMP  is  piecewise 
constant  and  right-continuous  with  left-hand  limits  everywhere.  Figure  4.1  shows 
a  sample  path  of  a  SMP,  along  with  its  embedded  Markov  renewal  sequence.  The 
process  starts  at  time  epoch  So  in  state  Y0,  where  So  =  0,  and  is  observed  again  at 
time  epoch  Si  to  be  in  state  Yx .  The  process  continues  to  evolve  and  is  observed  at 
discrete  time  epochs. 

Jf.,2  Instantaneous  Availability 

In  order  for  the  instantaneous  availability  to  be  calculated  for  a  satellite  degra¬ 
dation  process  modeled  as  a  SMP,  the  sojourn  time  distributions  for  each  state  of 
the  sample  space  S  must  be  known  as  well  as  the  state  transition  probabilities. 
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Figure  4.1  Sample  path  of  a  SMP. 

For  this  chapter,  the  probability  mass  function  a(t)  is  defined  as  it  was  for 
Chapter  3.  Recall  that  atj{t)  denotes  the  probability  that  the  process  is  in  state  j  at 
time  t,  defined  as 

aj(t)  =  p{x(t)=j}>  J  e  S’ 

and  the  probability  mass  function  ot(t)  is  a  row  vector  of  dimension  \S\,  where 

a(t)=  «i  (t)  oi2{t)  ...  a|S|(t)  • 

The  row  vector  a(0)  is  the  initial  probability  distribution  which  corresponds  to  the 
probabilities  associated  with  the  process  being  in  each  of  the  possible  states  in  S  at 
time  t  =  0,  where 

«(0)  =  a1(0)  a2(0)  ...  a|s|(0)  • 

The  matrix  G  (t),  consisting  of  elements 

Gijit)  =  P{ Xl  =  j,  S \  <  t\X0  =  i}  (4.3) 

defined  by  Eqn.  (4.2),  is  referred  to  as  the  semi-Markov  kernel.  A  SMP  is  completely 
described  by  its  initial  probability  distribution  a(0)  and  its  semi-Markov  kernel  G  (t), 
[30]. 

The  semi-Markov  kernel  is  directly  related  to  the  sojourn  time  distributions 
of  the  SMP.  Let  Tl3  denote  the  random  amount  of  time  the  process  stays  in  state  i 
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given  that  it  next  transitions  to  state  j.  The  random  variable  Ttl  has  a  c.d.f. 
where 

Hij(t)  =  P{Si  <  t\X0  =  i,Xi  =  j}. 

Let  Tt  denote  the  random  (unconditional)  amount  of  time  the  process  is  in  state  i. 
We  seek  the  c.d.f  H^t)  for  the  random  variable  T%.  By  conditioning  on  the  possible 
states  to  which  the  process  may  transition,  not  including  state  i, 


H,(t)  =  V  HiJ(t)P{Xl  =  j} 

j£S 

=  ^  t\Xo  =  b  =  j}P{X  1  =  j} 

jes 

=  Y,p{xi=j>S1<t\X0  =  i}. 

ieS 


Therefore, 

=  (4.4) 

jcs 

which  indicates  that  the  sojourn  time  distributions  are  composed  of  the  row  sums  of 
the  semi-Markov  kernel. 

The  semi-Markov  kernel  can  be  obtained  using  the  conditional  state  sojourn 
time  distributions  and  the  transition  probabilities  as 


Gij(t)  =  P{X1=j,S1<t \X0  =  i} 

=  P{ X1  =  j\X0  =  i}P{S1  <  t\X0  =  i,  X!  =  j} 

=  PijHijit).  (4.5) 

Solving  for  the  semi- Markov  kernel  matrix  in  this  manner  requires  knowledge  of  the 
c.d.f.  associated  with  the  sojourn  time  in  state  i  given  that  the  process  will  transition 
to  state  j.  This  also  requires  the  transition  probability  associated  with  the  process 
transitioning  from  state  %  to  state  j. 
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The  data  required  to  apply  this  methodology  to  a  satellite’s  degradation  process 
is  difficult  to  obtain  in  practice  and  as  such,  we  make  the  assumption  that  the  sojourn 
time  of  a  given  state  is  independent  of  the  subsequent  state  so  that 

Hi{t)  =  Hi:j(t),  for  all  j  ^  i. 

With  this  assumption,  we  may  write 


Gijit)  =  PijHiit).  (4.6) 

However,  as  more  historical  data  is  collected  regarding  the  distributions,  the  model 
may  be  updated  to  reflect  the  sojourn  time  distributions  for  specific  state  transitions. 
One  possible  way  to  incorporate  new  distribution  information  is  as  follows.  Let  151  —  1 
denote  the  total  number  of  states  to  which  the  process  can  transition,  not  including 
the  current  state.  In  general,  the  state  sojourn  time  distribution  is  given  as 


Hi{t)  =  P{Tt  <  t} 

=  1  -  P{Tt  >  t} 

=  1  -  P  {min(Tj;i,  Ti>2,  ■■■,  7i,|s|-i)  >  t} 

=  1  —  P {Tt.\  >  t,  Tt.2  >  t, . . . ,  Tit\S\-i  >  t}.  (4-7) 

Since  the  conditional  sojourn  time  distributions  are  independent, 

fl1(f)=i-nslj(i),  (4.8) 

where  Hij(t)  is  the  complement  of  Hi  jit). 

In  the  case  of  the  CTMC  model,  all  of  the  sojourn  time  distributions  are  expo¬ 
nential;  Hi{t)  is  simply  the  distribution  of  the  minimum  of  exponentially  distributed 
random  variables  which  is  exponential  with  rate  parameter  equal  to  the  sum  of  the 
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rates  associated  with  the  random  variables  T^i,  T^,  ■  •  •  ,T^s\-i-  The  rates  for  the 
state  holding  time  distributions  for  a  CTMC  correspond  to  the  diagonal  elements  of 
the  rate  matrix  Q,  where  the  negative  of  the  ith  diagonal  element  corresponds  to 
the  rate  of  the  ith  state  holding  time  distribution.  The  resulting  semi-Markov  kernel 
matrix  has  elements 


Gijit) 


0,  i  —  j 


(4.9) 


In  the  case  of  the  SMP  model,  the  state  holding  times  may  not  be  distributed 
exponentially;  therefore  the  transition  probabilities  plJ  may  not  be  known.  One  way 
to  approximate  the  semi-Markov  kernel  is  to  approximate  the  transition  probabilities 
by  assuming  that  a  statistical  sample  of  the  transition  rates  can  estimate  the  tran¬ 
sition  probability  as  if  the  distributions  of  the  state  holding  times  are  exponential. 
With  estimates  for  the  transition  rates,  we  can  approximate  the  semi-Markov  kernel, 
solve  a  system  of  equations  for  P  (t),  and  find  the  instantaneous  availability  by  using 
Eqn.  (3.18). 


Let  rij  denote  a  statistical  estimator  for  the  rate  of  transition  from  state  i  to 


state  j  given  by 


_N» (T) 

'  ij  T  ' 


(4.10) 


where  T  is  the  observation  period  and  iVy(T)  is  the  random  number  of  transitions 
from  state  i  to  state  j  in  time  period  T.  Note  that  the  observed  number  of  transitions 
and  the  period  of  observation  need  to  be  sufficiently  large  to  provide  a  reasonable 
statistical  estimator  for  the  true  rate  of  transition  [39].  Also,  we  may  take  as  the 
approximate  rate  of  transition  out  of  state  i 


\  ^ 

nr*  ry*  .  .  \  nr*  .  . 

1  %  '  ll  /  J  1  IJ  * 

j/* 


(4.11) 
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From  Eqn.  (4.10)  and  Eqn.  (4.11),  we  can  construct  a  matrix  R  which  represents  an 
estimator  for  the  rate  matrix  Q  where 

nj  «  qij  and  n  «  qt  . 


The  semi-Markov  kernel  can  now  be  approximated  by  applying  the  above  statistics 
to  Eqn.  (4.6)  where 


GiAt) 


i  A  j 


0, 


i  =  j. 


(4.12) 


Let  the  matrix  D(f)  be  a  diagonal  matrix  with  elements 


Hi(t),  i—j 

0,  i  A  j 


(4.13) 


The  following  equation  characterizes  the  relationship  between  P  (t),  G  (t),  and  D(f) 
[30]: 

P{t)  =  B{t)  +  G*P{t),  (4.14) 

where 

G*P(t)  —  [  G(t-u)dP(u)  (4.15) 

Jo 

is  the  convolution  of  G(f)  with  P(f). 

Let  f{t)  be  a  function  on  the  positive  real  line  and  let  f(s)  denote  the  Laplace- 
Sticltjes  transform  (LST)  of  f(t)  where  s  is  a  complex  variable.  The  Laplace-Stieltjes 
transform  is  defined  as 


f(s )  —  I  e  stdf(t),  t  >  0,  s  G  C,  3?(s)  >  0.  (4.16) 

Jo 

The  Laplace-Stieltjes  transform  of  matrix  P(f)  will  be  denoted  P(s)  where  its  (i,  j)th 
element  is  Pij(s). 
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An  important  theorem  relative  to  Laplace-Stieltjes  transforms  and  the  convo¬ 
lution  of  two  functions  states  that  given  two  functions  F(t)  and  G(t), 

LST{F*G(t)}  =  F(s)G(s). 

By  taking  the  LST  of  both  sides  of  Eqn.  (4.14), 

P(s)  =  D(s)  +  G(s)P(s).  (4.17) 


After  applying  matrix  algebra  to  Eqn.  (4.17),  the  Laplace  Sticltjes  transform  of  the 
probability  transition  matrix  P(s)  is  given  as 

P(s)  =  (I-G(S))-1D(s).  (4.18) 


By  properties  of  Laplace  and  Laplace-Stieltjes  transforms  [30], 

P *(s)  =  -P (s)  =  -(I  -  G(s))'1D(s).  (4.19) 

s  s 


Using  a  numerical  approximation  to  the  inverse  Laplace  transform  yields  a  solution 
for  the  probability  transition  matrix  P(t)  in  the  time  domain,  where 


P  (t)  =  C -1 


G(s))-  D(s) 


(4.20) 


Using  Eqn.  (4.20)  in  conjunction  with  Eqns.  (3.3)  and  (3.18),  the  satellite’s  instan¬ 
taneous  availability  may  be  obtained. 

The  LST  of  G (t)  does  not  exist  in  closed  form  for  sojourn  time  distributions 
that  are  power-tailed  and  for  some  that  are  heavy-tailed  [33].  A  sojourn  time  distri¬ 
bution  H-jit)  is  heavy-tailed  if 


lim 

f— XX) 


Hj(t  +  7) 
Hi(t) 


1,  7  >  0 


(4.21) 
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and  power-tailed  if 


(4.22) 


lim  Hi(t)  =  t  “,  a  >  0,  t  >  0. 

t— XX) 

For  these  types  of  c.d.f.s,  a  transform  approximation  method  can  be  applied  to 
approximate  G(s)  for  Eqn.  (4.20). 

The  following  sections  will  show  solution  methods  for  obtaining  the  Laplace 
transform  of  the  transition  probability  matrix  P  (t)  for  two-state  SMPs.  The  solution 
of  P*(s)  will  be  shown  for  two  cases  via  calculation  of  G(s)  and  D(s).  Each  case 
will  have  a  unique  distribution  for  its  state  holding  time.  An  approximation  method 
from  [18]  is  also  presented  for  the  Laplace  transform  of  Weibull  distributed  state 
holding  times,  as  the  transform  does  not  exist  in  closed  form. 

4-2.1  Example  1 

Let  the  state  holding  time  distributions,  mean  state  holding  estimates,  and 
the  associated  estimated  rates  for  each  state  of  S  be  given  as  in  Table  4.1.  Let  the 

Table  4.1  Example  1  holding  time  distributions  for  a  two-state  SMP. 


State 

Distribution 

Mean 

Ti 

1 

Gamma(fc,  A) 

k/X 

X/k 

2 

Uniform(a,  b ) 

(a  +  b) /2 

2/ (a  +  b ) 

transition  rate  approximations  from  state  i  to  state  j  be  given  as  rir 

To  calculate  P*(s),  solutions  for  G(s)  and  D(s)  are  substituted  into  Eqn.  (4.14) 
to  obtain  P(s).  Then  by  properties  of  Laplace  and  Laplace-Sticltjes  transforms, 
division  by  s  yields  P*(s).  First, 

Gu(t)  =  G22  (t)  =  0, 
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by  definition  of  the  semi-Markov  kernel,  which  implies  that 

Gfil(s)  =  G22(s)  =  0. 

For  a  two-state  SMP, 

rJi  =  y 
Ti 

since  the  only  transition  from  state  i  is  to  state  j  and  vice  versa.  Therefore, 

Gl2(t)  ~  —H^t)  =  Hi(t),  and  G2l(t)  «  =  H2(t), 

r  i  r2 

which  implies  that 

Gi2(s)  =  Hi(s),  and  G2i(s)  =  H2(s). 

By  definition  of  the  matrix  D  (t), 

D\2(t)  =  D2i{t)  =  0,  Duit)  =  1  —  Hi(t),  and  D22(t )  =  1  —  H2{t). 
The  Laplace-Stieltjes  transform  has  no  effect  on  a  constant.  Therefore, 
-Dii(s)  =  1  -  Hi(s),  and  D22(s)  =  1  —  H2{s). 


From  Eqn.  (4.16)  and  from  Table  4.1, 

From  these  results,  we  obtain  the  matrices  G(s)  and  D(s)  as 

0  (^r r ' 

G(s)  =  Kx+sJ 

-n^(e~sa  -  e~sb)  0 

s(b—a)  \  ' 
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and 

DM 

Therefore, 

P*M  =  i(l-GM)’1DM 

(e— -  e-") 
s(b  —  a) 

Numerical  methods  can  be  used  to  perform  the  inverse  Laplace  transform  to  obtain 
P(f)  in  the  time  domain. 

4-2.2  Example  2 

Let  the  state  holding  time  distributions,  mean  state  holding  estimates,  and 
the  associated  estimated  rates  for  each  state  of  S  be  given  as  in  Table  4.2.  From 


Table  4.2  Example  2  holding  time  distributions  for  a  two-state  SMP. 


State 

Distribution 

Mean 

n 

1 

Triangular(a,  m,  b) 

(a  +  m  +  b)/ 3 

3/(a  +  m  +  b) 

2 

Weibull(a,  A) 

r((a  +  l)/a)/A  A/r((a  +  l)/a) 

!-(£)• 

(1  _  (  A 

a  \A+s/  '  s(b— a) 


(e~sa—e~sb) " 
s(b—a) 

_  (e~SQ-e-sb) 
s(b—a) 


(AT 


->—sb\ 


s(b—a) 


Example  1,  it  was  shown  that  for  a  two-state  SMP,  P*(s)  is  simply  a  function  in 
terms  of  Hi(s)  and  H-zis).  From  Eqn.  (4.16)  and  Table  4.2, 

^  2  /  g— sa  g—sm  g— sb  ^—sm  \ 

s2  \ v(5  —  a)(m  —  a)  (b  —  a)  (b  —  m) )  ' 

A  closed  form  solution  does  not  exist  for  H2(s)  due  to  the  fact  that  the  Weibull 
distribution  is  a  heavy-tailed  distribution  and  does  not  possess  an  analytical  Laplace 
transform.  Therefore,  a  numerical  approximation  to  the  Laplace  transform  must  be 
used  to  obtain  P*(s). 
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The  Transform  Approximation  Method  (TAM)  was  developed  for  the  purpose 
of  finding  an  approximation  to  heavy-tailed  distributions  and  is  explained  in  numer¬ 
ous  articles  including  [33]  and  [18].  A  discrete  approximation  is  performed  on  N 
discrete  points  in  the  domain  of  some  c.cl.f.  F(x),  where  each  point  is  referred  to 
as  a  TAM  sample  and  is  denoted  x(i),  i  =  1,  2, . . . ,  N.  Each  of  these  points  has  an 
equal  probability  mass,  resulting  in  a  cumulative  distribution  estimate 

H^))  = 

Then  using  a  method  from  [22], 

1  N 

2=1 

From  [18],  the  TAM  samples  for  the  Weibull  distribution  are  given  by 

x(i)  =  0  [“  1"  (l  VTl)]  ”  ’ 

resulting  in  the  approximation 

n  r  f  ,  x-,i' 

H2(s)  «  sH*(s)  =  ^J^exp  ~SP  _ln(1_ArTl) 

v  i=l  L  L  V  /  J  _ 

From  Hi(s)  and  the  approximation  for  ^(s),  P*(s)  can  be  obtained  through 
algebra.  Numerical  methods  can  then  be  used  to  perform  the  inverse  Laplace  trans¬ 
form  to  obtain  P  (t)  in  the  time  domain. 

4-3  Limiting  Availability 

The  methods  for  obtaining  the  limiting  availability  for  a  satellite’s  degradation 
process  modeled  as  a  SMP  are  given  below  for  comparison  purposes.  Let  A(k^  and 
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A^k\t)  denote  the  limiting  availability  and  the  instantaneous  availability  for  satellite 
k  at  time  t,  respectively.  Similar  to  Chapter  3,  the  limiting  availability  is  solved  by 
first  obtaining  the  probability  distribution  for  the  degradation  process  in  its  asymp¬ 
totic  regime,  then  by  summing  up  the  long-run  probabilities  associated  with  states 
of  availability. 

Let  p-C  denote  a  row  vector  consisting  of  probabilities  }  defined  as 
limP(XW(i)  =  at),  tneS,  i  =  l,2,...,\S\. 

t—*oo 

The  probability  represents  the  long-run  probability  the  satellite  is  found  in  state 
i.  For  a  SMP  that  is  irreducible  and  recurrent,  the  transition  probability  matrix 
P  =  G  (oo).  Let  tt^  denote  a  row  vector  of  probability  values,  where 

n(k)  _  7r(fc)G(o0). 


The  system  of  equations 


7r(fc)P(fc)  =  0 
7 r^e  =  1 

is  solved  for  nl'k\  where  P(fc)  is  the  transition  probability  matrix  associated  with 
satellite  k ,  0  denotes  a  row  vector  of  zeros  and  e  denotes  a  row  vector  of  ones.  As  was 
stated  earlier  in  the  chapter,  since  the  transition  probabilities  for  the  SMP  satellite 
model  are  not  known,  we  use  statistics  to  approximate  them.  From  Eqn.  (4.10)  and 
Eqn.  (4.11),  an  approximation  for  7r®  can  be  solved  from  the  system  of  equations 

=  o 

7r^e  =  1 
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where  the  matrix  RjC  represents  the  statistical  estimator  for  the  rate  matrix  of 
satellite  k. 


After  solving  for  7r^k\  the  probabilities  of  row  vector  can  be  solved  by  the 
equation 

pf]  =  lirn  P(X(fc)(t)  =  i)  = 


(k) 

K  ki 


t— >oo 


esav 

where  /q  is  the  expected  state  holding  time  of  state  i.  Finally,  the  limiting  availability 


for  satellite  k  is  obtained  by  summing  the  long-term  probabilities  for  the  states  of 
availability  resulting  in 


\s'\ 

A{k)  =  y 


i= 1 


4-4  Model  Summary 

This  chapter  has  presented  a  feasible  extension  to  the  formulation  of  the  con¬ 
stellation  functional  capability  measure.  The  satellite  instantaneous  availability  mea¬ 
sure  was  generalized  by  modeling  the  degradation  process  as  a  semi-Markov  process. 
A  methodology  was  described  for  acquiring  the  state  probabilities  at  a  given  time 
when  the  holding  time  distributions  are  generally  distributed.  A  brief  overview  of 
SMPs  was  provided  along  with  examples  of  computing  the  Laplace  transform  of  the 
probability  transition  matrix  for  a  SMP  with  distinct  sojourn  time  distributions. 
Finally,  a  methodology  was  described  for  obtaining  the  limiting  availability  of  a 
satellite. 

In  the  next  chapter,  we  illustrate  the  methodologies  presented  in  Chapters 
3  and  4  with  three  numerical  examples.  Each  example  consists  of  calculating  the 
constellation  functional  capability  measure  via  the  satellite  instantaneous  availability 
measures  and  value  scores.  Results  from  these  examples  are  compared  with  those 
obtained  via  Monte  Carlo  simulation. 
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5.  Numerical  Results 


This  chapter  presents  a  few  examples  to  illustrate  the  methodology  developed 
in  Chapters  3  and  4.  For  each  example,  the  sample  space  is  described  for  the  satellites 
contained  within  their  parent  constellation,  then  the  instantaneous  availability  is 
calculated  for  each  satellite  based  on  the  distribution  of  the  failure  and  repair  rates; 
the  limiting  availability  is  also  calculated  in  the  third  example.  Finally,  each  satellite 
is  given  a  value  score  based  on  a  set  of  defined  attributes  relating  its  contribution  to 
the  constellation’s  mission. 

For  these  examples,  the  instantaneous  availability  and  value  scores  are  com¬ 
puted  using  notional  data;  however,  the  examples  are  carefully  constructed  to  repre¬ 
sent  satellite  design  lifetimes  from  open  literature  resources.  The  results  will  demon¬ 
strate  the  benefits  associated  with  the  computational  speed  of  using  analytical  meth¬ 
ods  rather  than  simulation  models  to  calculate  satellite  availability  as  well  as  the 
errors  prone  to  measuring  a  satellite’s  limiting  availability  over  its  instantaneous 
availability. 

5.1  Description  of  Experiments 

In  this  section  we  describe  the  means  by  which  we  compute  the  functional 
capability  measures  for  each  of  the  given  examples.  The  comparison  between  the 
analytical  and  simulation  methods  is  based  on  the  computation  of  the  satellite’s 
instantaneous  availability.  The  instantaneous  availability  measures  for  the  analytical 
results  were  calculated  using  the  CTMC  methodology  of  Chapter  3  for  example  1, 
and  the  SMP  methodology  of  Chapter  4  for  the  remaining  examples.  Both  the 
analytical  results  and  the  simulations  were  performed  using  MatLab  software  on  a 
Dell  Precision  670  computer  with  dual  Intel  Xeon  2.8  GHz  processors  with  3  GB 
of  RAM.  The  satellite  value  scores  and  the  final  constellation  functional  capability 
measure  were  computed  using  Eqns.  (3.28)  and  (3.1),  respectively,  of  Chapter  3. 


5-1 


The  first  steps  of  computing  the  analytical  instantaneous  availability  focused 
on  constructing  the  rate  matrix  Q,  where  a  subprogram  was  written  for  this  specific 
purpose.  The  subprogram  first  produces  all  of  the  state  space  elements,  then  de¬ 
termines  which  elements  correspond  to  transitions  of  function  failures  and  repairs. 
For  the  CTMC  case,  the  rates  were  obtained  from  the  failure  and  repair  rates  of  the 
functions.  For  the  SMP  case,  the  distribution  of  the  state  holding  times  were  used  to 
construct  the  rates.  By  using  the  negative  reciprocals  of  the  expected  state  holding 
times  for  the  diagonal  entries,  the  remaining  entries  were  constructed  to  make  the 
rows  sum  to  zero. 

Next,  the  Laplace  transform  of  the  transition  probability  matrix  P(f)  was 
computed  from  Eqn.  (4.19).  Analytical  expressions  of  the  Laplace  transforms  of  the 
distribution  functions  were  used  to  calculate  the  transition  probability  matrix  from 
the  matrices  G(s)  and  D(s).  In  the  case  where  an  analytical  expression  did  not 
exist  for  the  distribution,  an  approximation  to  the  Laplace-Stieltjes  transform  was 
performed  using  the  transform  approximation  method  discussed  in  Chapter  4  with 
N  =  1000.  Numerical  Laplace  transform  inversion  was  then  used  to  obtain  a  time 
domain  solution  for  the  transition  probability  matrix.  Each  inverse  Laplace  trans¬ 
form  was  computed  using  the  algorithm  of  Abate  and  Whitt  [1] .  Once  the  transition 
probability  matrix  was  obtained,  Eqn.  (3.3)  was  solved  to  obtain  the  probability  mass 
function,  and  Eqn.  (3.18)  was  solved  to  obtain  the  satellite  instantaneous  availability. 

A  simulation  model  was  also  developed  and  executed  to  acquire  the  instan¬ 
taneous  availability  results  for  each  satellite  in  the  constellation.  These  simulation 
results  were  compared  to  the  analytical  results  in  terms  of  processing  time,  as  the 
actual  numbers  were  indistinguishable  in  most  cases.  The  simulation  consisted  of 
running  multiple  replications  of  the  degradation  process  evolving  randomly  up  to 
time  t,  the  time  at  which  the  instantaneous  availability  is  evaluated.  The  process 
began  in  state  1  and  remained  in  this  state  for  a  random  amount  of  time  determined 
by  a  random  number  generator  corresponding  to  the  state  holding  distribution  of 
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state  1.  A  subfunction  randomly  transitions  the  process  to  its  next  state  by  incor¬ 
porating  the  rate  matrix  Q.  The  process  remains  in  the  next  state  for  a  random 
amount  of  time  determined  by  a  random  number  generator  corresponding  the  that 
state’s  holding  time  distribution. 

Transitions  from  one  state  to  another  were  simulated  with  the  same  rate  matrix 
used  in  the  analytical  computations.  After  each  of  the  0.5  million  replications,  a 
counter  records  the  number  of  times  the  process  is  found  in  each  of  the  possible 
states  at  time  t.  A  probability  mass  function  is  constructed  using  these  tallies.  From 
this  probability  mass  function,  the  states  of  availability  are  summed  to  obtain  the 
simulated  instantaneous  availability  of  the  satellite  at  time  t. 

The  selected  value  scores  for  the  three  examples  were  based  on  the  method¬ 
ology  presented  in  Chapter  3.  Each  example  uses  the  same  set  of  attributes  to 
score  the  satellites.  Values  for  constructing  the  attribute  value  functions,  weights, 
and  ultimately  the  satellite  score,  were  notional.  Using  Eqn.  (3.1),  the  constella¬ 
tion  functional  capability  measure  was  computed  with  each  satellite’s  instantaneous 
availability  and  value  score. 

5.2  Navstar  GPS  Constellation 

The  U.S.  Military  relies  heavily  on  the  Navstar  Global  Positioning  System 
(GPS)  to  provide  accurate  time,  location,  and  velocity  data  to  mobile  platforms 
such  as  aircraft  and  munitions  [38].  The  performance  of  this  constellation  depends 
on  the  capability  of  its  24  satellites  to  perform  mission  essential  navigation  functions. 
By  calculating  each  satellite’s  instantaneous  availability  and  scoring  each  satellite 
based  on  its  operational  contribution  to  the  GPS  constellation  mission,  the  functional 
capability  of  this  constellation  will  be  assessed  at  t  —  9  years. 

Suppose  that  the  constellation  consists  of  16  Block  IIA  satellites,  and  8  Block 
HR  satellites,  where  the  Block  HR  satellites  are  newer,  are  designed  with  the  latest 
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technology,  and  are  built  to  last  longer  than  the  Block  IIA  satellites.  Each  satellite 
is  independent  of  all  other  satellites  and  possesses  three  independent  functions  that 
have  exponential  failure  and  repair  rates,  the  satellite’s  clock  being  the  first  function, 
the  computer  being  the  second,  and  the  transceiver  the  third.  A  subject  matter 
expert  familiar  with  the  constellation’s  mission  and  how  its  satellites  contribute  to 
the  mission  is  used  to  determine  a  subset  of  the  state  space  indicative  of  the  satellite 
being  available.  For  demonstrative  purposes,  suppose  that  in  order  for  these  satellites 
to  be  considered  available,  all  three  functions  must  be  available. 

5.2.1  Instantaneous  Availability 

The  state  space  for  each  satellite  will  be  identical  and  are  defined  in  Table  5.1. 
Recall  from  Eqn.  (3.2),  the  random  variable  indicating  the  state  of  a  satellite  is  a 
vector  of  elements  where  a  value  of  1  refers  to  the  function  being  available  and  a  value 
of  0  refers  to  the  function  being  unavailable.  Figure  5.1  shows  the  transition  rate 

Table  5.1  Satellite  states  for  the  Navstar  constellation  example. 


State 

Combination 

State 

Combination 

1 

(1,1,1) 

5 

(0,0,1) 

2 

(o,i,i) 

6 

(0,1,0) 

3 

(1,0,1) 

7 

(1,0,0) 

4 

(1,1,0) 

8 

(0,0,0) 

diagram  for  each  satellite  in  this  example,  showing  the  possible  transitions  that  can 
occur  between  states.  The  satellites  are  listed  in  Table  5.2  with  their  corresponding 
initial  probability  distributions  and  their  years  of  initialization. 

In  this  example,  at  the  initial  time  satellites  17  through  20  start  their  useful 
life  in  the  constellation,  their  status  cannot  be  determined  with  certainty;  there  is  a 
small  probability  associated  with  hireling  the  satellite  in  state  3,  corresponding  to  a 
0.02  probability  that  the  computer  is  not  working  at  time  to-  This  example  illustrates 
how  the  model  can  incorporate  a  situation  where  the  initial  status  of  the  satellite 
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Figure  5.1  Transition  rate  diagram  for  the  Navstar  constellation  ex¬ 
ample. 

functions  cannot  be  determined  with  certainty.  Regardless  of  the  interpretation  of 
the  satellite’s  initialization  time  to,  the  initial  probability  distribution  contains  the 
probabilities  associated  with  the  satellite  being  in  a  given  state  at  to-  If  a  satellite 
check-out  can  be  performed  at  the  time  a  satellite  starts  its  life  in  the  constellation, 
then  the  state  of  the  satellite  will  be  known  with  certainty. 

Table  5.3  lists  each  function’s  failure  and  repair  distribution.  Notice  that  the 
computer  is  the  only  function  considered  to  be  repairable  in  this  example. 


Table  5.2  Satellite  information  for  the  Navstar  constellation  example. 


Satellites 

Type 

Initial  Probability  Distribution 

Time  Initialized  (t0) 

1-4 

Block  IIA 

(1,0, 0,0, 0,0, 0,0) 

0 

5-16 

Block  IIA 

(1,0, 0,0, 0,0, 0,0) 

2 

17-20 

Block  HR 

(0.98,0,0.02,0,0,0,0,0) 

4 

21-24 

Block  HR 

(1,0, 0,0, 0,0, 0,0) 

6 

A  graph  of  the  instantaneous  availability  measures  is  given  for  each  of  the 
four  groups  of  satellites  in  Figure  5.2,  calculated  using  the  methods  of  Chapter  3. 
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Tabic  5.3  Function  information  for  the  Navstar  constellation  example. 


Satellites 

Function 

Failure  Distribution 

Repair  Distribution 

1-4 

Clock 

Computer 

Transceiver 

Exp(0.0291) 

Exp(0.0503) 

Exp(0.0230) 

0 

Exp(3.3300) 

0 

5-16 

Clock 

Exp(0.0121) 

0 

Computer 

Exp(0.0421) 

Exp(4.8216) 

Transceiver 

Exp(0.0189) 

0 

17-20 

Clock 

Exp(0.0121) 

0 

Computer 

Exp(0.0303) 

Exp(4.9216) 

Transceiver 

Exp(0.0143) 

0 

21-24 

Clock 

Exp(0.0103) 

0 

Computer 

Exp(0.0282) 

Exp(4.9216) 

Transceiver 

Exp(0.0118) 

0 

The  simulated  availability  measures  were  not  included  in  this  figure  as  they  were 
indistinguishable  from  the  analytical  measures.  Notice  that  satellites  17  through  20 
start  with  0.98  availability  at  their  initial  times,  corresponding  to  the  uncertainty 
associated  with  the  state  these  satellites  may  be  in  at  time  t0.  Table  5.4  summarizes 
the  instantaneous  availability  measures  for  each  of  the  satellites  at  t  —  9  years. 
These  measures  take  into  account  the  amount  of  time  the  satellite  has  been  in  orbit 
by  calculating  its  instantaneous  availability  based  on  its  time  of  initialization,  t0.  The 
average  time  for  the  analytical  availability  computation  was  0.01  minutes  versus  2.62 
minutes  for  the  simulation  method.  Thus,  the  computation  increased  by  a  factor  of 
over  243  for  the  simulation  results.  The  processing  time  for  the  simulation  may 
be  negligible  in  this  case,  however,  the  analytical  solution  may  be  a  more  feasible 
option  for  larger  constellations  containing  satellites  with  more  than  three  functions. 
Also,  to  obtain  a  simulation  result  with  statistical  confidence  would  require  more 
replications,  increasing  the  computational  effort. 

The  limiting  availability  for  all  of  the  24  satellites  is  equal  to  zero  as  some 
of  the  required  functions  cannot  be  repaired  when  they  fail.  When  the  satellite 
experiences  a  failure  with  its  clock  or  with  its  transceiver,  it  is  considered  unavailable 
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Figure  5.2  Comparison  of  instantaneous  availability  for  Navstar  satellites 
computed  analytically. 

and  stays  in  this  state  for  the  remainder  of  its  lifetime.  Since  the  limiting  availability 
is  the  proportion  of  uptime  over  total  time,  the  numerator  stays  constant  while  the 
denominator  approaches  infinity,  resulting  in  a  limiting  availability  equal  to  zero. 

5.2.2  Value  Scores 

Each  satellite  is  assumed  to  be  judged  by  a  SME  based  on  three  independent  at¬ 
tributes:  (1)  contribution  to  mission,  (2)  level  of  technology,  and  (3)  signal  accuracy. 
Attributes  1  and  2  were  developed  by  [35].  These  attributes  are  non-specific  and 
could  be  applied  to  a  number  of  different  constellations.  The  normalized  attribute 
scales  are  given  in  Table  5.5. 

For  each  satellite,  a  value  score  is  calculated  using  the  methods  of  Chapter  3 
by  creating  a  value  function  and  weight  for  each  attribute,  then  having  a  SME  score 
each  attribute  of  each  satellite.  The  value  functions  and  attribute  weights  for  the 
three  attributes  are  given  in  Table  5.6  and  are  obtained  from  computing  the  risk 
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Table  5.4  Sample  availability  measures  for  the  Navstar  constellation  example  evaluated 
at  t  =  9  years. 


Satellites 

to 

t-t0 

A(t 

Analytical 

—  to) 

Simulated 

Difference 

1-4 

0 

9 

0.616380 

0.615982 

0.000398 

5-16 

2 

7 

0.797963 

0.797916 

0.000047 

17-20 

4 

5 

0.870979 

0.871156 

0.000177 

21-24 

6 

3 

0.930518 

0.930296 

0.000222 

Table  5.5  Normalized  attribute  scales  for  the  Navstar  constellation  example. 


Attribute 

Measure  (y) 

Description 

Contribution 

0.00 

No  Direct  Impact 

to  Mission 

0.25 

Secondary  Impact  to  Mission 

0.50 

Impacts  Mission 

0.75 

Strongly  Impacts  Mission 

1.00 

Critical  to  Mission 

Level  of 

0.00 

Primitive 

Technology 

0.25 

30  Year  Old  Technology 

0.50 

20  Year  Old  Technology 

0.75 

10  Year  Old  Technology 

1.00 

State  of  the  Art 

Signal 

0.00 

No  Usable  Signal 

Accuracy 

0.33 

Low  Accuracy 

0.66 

Medium  Accuracy 

1.00 

High  Accuracy 

attitude  constants  by  solving  Eqn.  (3.30)  with  Newton’s  method.  These  scores  are 
time  dependent,  as  the  attribute  weights  are  assigned  at  the  time  of  evaluation  based 
on  the  preference  of  the  attributes  to  the  constellation’s  mission  at  time  t.  Table  5.7 
summarizes  the  results  associated  with  the  value  scores  for  each  of  the  satellites.  As 
in  Chapter  3,  Vi{yi)  is  the  value  function  for  attribute  i  evaluated  at  score  yl  and 
V^(9)  is  the  score  for  satellite  k  at  time  t  —  9  years,  calculated  using  Eqn.  (3.28) 
with  the  satellite’s  value  functions  and  the  attribute  weights. 
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Table  5.6  Value  functions  and  attribute  weights  for  the  Navstar  constellation  example. 


Attribute  (i) 

y'i 

Q 

Vi{yi) 

Wi(9) 

1 

0.60 

0.822163 

1— exp(0.822163t/i) 
l—exp(0. 822163) 

0.31 

2 

0.65 

1.278652 

1— exp(1.278652j/i) 
l—exp(l. 278652) 

0.29 

3 

0.66 

1.376696 

1  —  exp(l  .376696?/; ) 
l—exp(l. 376696) 

0.40 

5. 2. 3  Functional  Capability 

Based  on  the  instantaneous  availability  measurements  along  with  the  value 
scores  for  each  satellite  in  the  constellation,  the  constellation  functional  capability 
measure  <3>(9)  is  calculated  to  be  0.710162  (0.710081  using  simulation  results). 


5.3  Milstar  Satellite  Communication  Constellation 

The  Milstar  satellite  communication  system  provides  essential  communication 
capabilities  to  the  U.S.  military  through  its  five-satellite  constellation  [37].  The 
performance  of  this  constellation  is  dependent  on  the  capability  of  its  5  satellites 
to  perform  missions  relating  to  secure,  global  communication.  By  calculating  each 
satellite’s  instantaneous  availability  and  scoring  each  satellite  based  on  its  opera¬ 
tional  contribution  to  the  Milstar  constellation  mission,  the  functional  capability  of 
this  constellation  will  be  assessed  at  t  —  9  years  using  the  methodology  described  in 
Chapters  3  and  4. 

For  this  example,  suppose  that  the  constellation  consists  of  5  identical  and 
independent  satellites.  Each  satellite  consists  of  an  independent  function  for  each 
of  the  four  types  of  information  it  processes:  voice,  data,  teletype,  and  facsimile. 
Further  suppose  that  in  order  for  these  satellites  to  be  considered  available,  the  Erst 
two  functions  must  be  available;  assume  that  if  the  teletype  and  facsimile  functions 
become  unavailable,  the  availability  of  the  data  function  keeps  the  satellite  in  an 
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Tabic  5.7  Value  score  results  for  the  Navstar  constellation  example. 


Satellite 

Vi 

V-2 

V3 

Mv  i) 

V2M) 

My 3) 

V(9) 

1 

0.84 

0.79 

0.98 

0.780095 

0.673654 

0.963672 

0.822658 

2 

0.86 

0.79 

0.94 

0.806027 

0.673654 

0.893949 

0.802808 

3 

0.84 

0.79 

0.86 

0.780095 

0.673654 

0.765511 

0.743393 

4 

0.84 

0.79 

0.95 

0.780095 

0.673654 

0.911022 

0.801598 

5 

0.96 

0.83 

0.97 

0.942283 

0.729252 

0.945880 

0.881943 

6 

0.86 

0.83 

0.98 

0.806027 

0.729252 

0.963672 

0.846820 

7 

0.95 

0.83 

0.94 

0.928148 

0.729252 

0.893949 

0.856789 

8 

0.97 

0.83 

0.99 

0.956535 

0.729252 

0.981711 

0.900693 

9 

0.98 

0.83 

0.99 

0.970904 

0.729252 

0.981711 

0.905148 

10 

0.94 

0.83 

0.98 

0.914128 

0.729252 

0.963672 

0.880332 

11 

0.86 

0.83 

0.99 

0.806027 

0.729252 

0.981711 

0.854036 

12 

0.95 

0.83 

0.89 

0.928148 

0.729252 

0.812026 

0.824019 

13 

0.83 

0.83 

0.96 

0.767288 

0.729252 

0.928331 

0.820675 

14 

0.94 

0.83 

0.98 

0.914128 

0.729252 

0.963672 

0.880332 

15 

0.86 

0.83 

0.99 

0.806027 

0.729252 

0.981711 

0.854036 

16 

0.95 

0.83 

0.96 

0.928148 

0.729252 

0.928331 

0.870541 

17 

0.83 

0.94 

0.99 

0.767288 

0.897656 

0.981711 

0.890864 

18 

0.99 

0.94 

0.99 

0.985392 

0.897656 

0.981711 

0.958476 

19 

0.89 

0.94 

0.98 

0.845734 

0.897656 

0.963672 

0.907967 

20 

0.96 

0.94 

0.97 

0.942283 

0.897656 

0.945880 

0.930780 

21 

0.98 

0.97 

0.99 

0.970904 

0.947847 

0.981711 

0.968540 

22 

0.97 

0.97 

0.98 

0.956535 

0.947847 

0.963672 

0.956870 

23 

0.99 

0.97 

0.99 

0.985392 

0.947847 

0.981711 

0.973032 

24 

0.98 

0.97 

0.99 

0.970904 

0.947847 

0.981711 

0.968540 

available  status.  Also,  suppose  that  the  only  function  that  can  be  repaired  on  the 
satellite  is  the  voice  function,  provided  that  the  data  function  is  still  available. 


5.3.1  Instantaneous  Availability 

The  state  space  for  each  satellite  will  be  identical  and  is  defined  in  Table 
5.8.  Notice  that  states  1,  4,  5,  and  11  are  states  of  availability  for  this  satellite, 
corresponding  to  having  both  the  voice  and  data  functions  available.  Figure  5.3 
shows  the  possible  transitions  that  can  occur  between  states.  The  satellites  are 
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Tabic  5.8  Satellite  states  for  the  Milstar  constellation  example. 


State 

Combination 

State 

Combination 

1 

(1,1, 1,1) 

9 

(1,0, 0,1) 

2 

(0,1, 1,1) 

10 

(1,0, 1,0) 

3 

(1,0, 1,1) 

11 

(1,1, 0,0) 

4 

(1,1, 0,1) 

12 

(0,0,0, 1) 

5 

(1,1, 1,0) 

13 

(0,0, 1,0) 

6 

(0,0, 1,1) 

14 

(0,1, 0,0) 

7 

(0,1, 0,1) 

15 

(1,0, 0,0) 

8 

(0,1, 1,0) 

16 

(0,0, 0,0) 

listed  in  Table  5.9  with  their  corresponding  initial  probability  distribution,  and  their 
year  of  initialization. 

Table  5.10  lists  the  state  holding  time  distributions  for  each  of  the  states  for 
this  example. 


Table  5.9  Satellite  information  for  the  Milstar  constellation  example. 


Satellite 

Initial  Probability  Distribution 

Time  Initialized  (to) 

1 

(1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0) 

0 

2 

(1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0) 

4 

3 

(1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0) 

6 

4 

(1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0) 

7 

5 

(1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0) 

8 
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Figure  5.3  Transition  rate  diagram  for  the  Milstar  constellation  ex¬ 
ample. 

Let  the  statistical  estimates  for  the  transition  rates  be  given  in  the  matrix  R. 
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0.000 
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0.000 
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0.000 

0.000 
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0.000 

0.000 
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0.000 
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0.000 
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0.000 

0.000 

0.000 
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Table  5.10  State  holding  time  distributions  for  the  Milstar  constellation  example. 


State 

Distribution 

State 

Distribution 

1 

Gamma(3.80,0.37) 

9 

Triangular  (1.80, 2. 10,2.30) 

2 

Gamma(2. 10,0.50) 

10 

Triangular  (1.61, 2. 20, 3. 20) 

3 

Gamma(2.20,0.41) 

11 

Triangular(1.10,1.83,2.40) 

4 

Gamma  ( 3 . 2 0 , 0 . 42 ) 

12 

Triangular(1.70,2.50,4.51) 

5 

Weibull(1.20,0.22) 

13 

Uniform(0.05,2.30) 

6 

Weibull(1.80,0.42) 

14 

Uniform(0.13,2.10) 

7 

Weibull(1.90,0.31) 

15 

Uniform(0.03,2.40) 

8 

Weibull(1.20,0.33) 

16 

Uniform(0.51,0.93) 

Since  the  voice  function  is  the  only  function  that  can  be  repaired,  the  approximation 
R  to  the  rate  matrix  Q  is  sparse  below  its  diagonal.  (A  satellite  consisting  of 
functions  that  cannot  be  repaired  results  in  a  lower  triangular  rate  matrix.) 

A  graph  of  the  instantaneous  availability  compared  to  the  simulated  availability 
for  the  satellites  in  this  constellation  is  given  in  Figure  5.4,  where  the  instantaneous 
availability  is  calculated  using  Eqns.  (4.6)-(4.20)  of  Chapter  4,  and  Eqns.  (3.3)  and 
(3.18)  of  Chapter  3.  Note  that  the  instantaneous  availability  measures  for  each  of 
the  5  satellites  can  be  derived  from  the  same  analysis  by  accounting  for  the  year 
of  initialization.  Table  5.11  summarizes  the  instantaneous  availability  measures  for 
each  of  the  satellites  at  t  —  9  years.  The  average  time  for  the  analytical  availability 
computation  was  1.16  minutes  versus  10.47  minutes  for  the  simulation  method.  Thus, 
the  computation  increased  by  a  factor  of  over  9  for  the  simulation  results. 


Table  5.11  Sample  availability  measures  for  the  Milstar  constellation  example  evaluated 
at  t  =  9  years. 


Satellite 

to 

t  - 10 

A(t 

Analytical 

—  to) 

Simulated 

Difference 

1 

0 

9 

0.722399 

0.740352 

0.017953 

2 

4 

5 

0.903839 

0.912136 

0.008297 

3 

6 

3 

0.961779 

0.965656 

0.003877 

4 

7 

2 

0.980183 

0.981590 

0.001407 

5 

8 

1 

0.991859 

0.992288 

0.000429 
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Figure  5.4  Comparison  of  instantaneous  availability  for  the  Milstar  satel¬ 
lites  computed  analytically  (solid  line)  and  simulated  availabil¬ 
ity  (dotted  line). 

As  was  the  case  with  the  Navstar  constellation  example,  the  limiting  availabil¬ 
ity  for  all  of  the  Milstar  constellation  satellites  is  equal  to  zero  as  there  exist  functions 
that  cannot  be  repaired  which  are  necessary  for  the  satellite  to  be  considered  avail¬ 
able.  For  example,  once  the  satellite  reaches  state  16,  it  will  never  leave  state  16  as 
there  are  no  repairs  that  will  take  place.  Since  the  data  function  is  required  for  the 
satellite  to  be  classified  as  available,  and  the  data  function  cannot  be  repaired,  the 
satellite’s  availability  is  limited  by  the  life  of  the  data  function. 

5.3.2  Value  Scores 

The  value  scores  of  the  Milstar  constellation  will  be  calculated  using  the  same 
set  of  attributes  that  were  used  to  score  the  satellites  in  the  Navstar  constellation 
example.  The  value  functions  and  attribute  weights  for  the  three  attributes  are  given 
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in  Table  5.19,  and  Table  5.20  summarizes  the  results  associated  with  the  value  scores 
for  each  of  the  satellites. 


Table  5.12  Value  functions  and  attribute  weights  for  the  Milstar  constellation  example. 


Attribute  (i) 

V'i 

Ci 

Vi{yi) 

Wi(  9) 

1 

0.73 

2.164803 

1  — exp  (2 . 1 64803j/i ) 
l-exp(2. 164803) 

0.57 

2 

0.64 

1.183190 

l-exp(1.183190j/i) 
l-exp(l. 183190) 

0.29 

3 

0.55 

0.402692 

1— exp(0.4026924/i) 
l—exp(0. 402692) 

0.14 

Table  5.13  Value  score  results  for  the  Milstar  constellation  example. 


Satellite 

IJ\ 

2/2 

2/3 

Mvi) 

^2(2/2) 

Vs  (2/3) 

V(9) 

1 

0.94 

0.73 

0.81 

0.862401 

0.605795 

0.777793 

0.776140 

2 

0.96 

0.83 

0.94 

0.906296 

0.737345 

0.927984 

0.860337 

3 

0.93 

0.90 

0.86 

0.841156 

0.839142 

0.834630 

0.839658 

4 

0.98 

0.94 

0.95 

0.952134 

0.901211 

0.939866 

0.935649 

5 

0.96 

0.96 

0.97 

0.906296 

0.933364 

0.963775 

0.922193 

5.3.3  Functional  Capability 

Based  on  the  instantaneous  availability  measurements  along  with  the  value 
scores  for  each  satellite  in  the  constellation,  the  constellation  functional  capability 
measure  $(9)  is  calculated  to  be  0.795529  (0.800737  using  simulation  results). 

5-4  Defense  Meteorological  Satellite  Constellation 

The  U.S.  Air  Force  collects  weather  data  from  a  constellation  of  two  mete¬ 
orological  satellites  in  polar  orbits  [36].  These  satellites  have  multiple  functions 
associated  with  performing  their  mission.  The  satellite’s  primary  system  collects  vi¬ 
sual  and  infrared  imagery  data.  There  are  also  secondary  sensors  on  the  satellite,  one 
which  “measure [s]  atmospheric  vertical  profiles  of  moisture  and  temperature”  and 
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another  which  measures  the  ionosphere  to  determine  its  impact  on  other  space-based 
systems  [36:1].  We  will  consider  the  satellites  in  this  constellation  to  be  independent 
and  satellite  1  has  four  functions  defined  in  Table  5.14  and  satellite  2  has  only  the 
first  three  of  these  functions.  Each  satellite  will  be  considered  available  if  all  of  its 
functions  are  available.  Also,  all  of  the  functions  will  be  considered  repairable.  The 
instantaneous  availability  and  value  score  will  be  calculated  for  both  satellites  of  this 
constellation  to  compute  the  constellation  functional  capability  measure  at  t  —  10 
years. 


Table  5.14  Satellite  functions  for  the  meteorological  constellation  example. 


Function 

Description 

1 

Visual  imagery  sensor 

2 

Infrared  imagery  sensor 

3 

Vertical  moisture  and  temperature  profile  sensor 

4 

Ionosphere  sensor 

5-4-1  Instantaneous  Availability 

The  state  space  for  satellite  1  is  the  same  as  the  state  space  for  the  satellites  in 
the  Milstar  constellation  example  (Table  5.8).  The  state  space  for  satellite  2  is  the 
same  as  the  state  space  for  the  satellites  in  the  Navstar  constellation  example  (Table 
5.1).  The  transition  rate  diagram  for  satellite  1  is  similar  to  the  diagram  in  Figure 
5.3  for  the  Milstar  constellation  example.  The  transition  rate  diagram  for  satellite  2 
is  similar  to  the  diagram  in  Figure  5.1  for  the  Navstar  constellation  example.  Since 
all  of  the  functions  can  be  repaired  on  the  meteorological  constellations  satellites,  all 
of  the  states  communicate,  resulting  in  transition  rate  diagrams  with  two-direction 
arcs. 

Both  of  these  satellites  are  assumed  to  have  initial  probability  distributions 
where  the  probability  of  being  in  state  1  is  equal  to  1.0.  Both  satellites  will  also 
be  assumed  to  have  initialized  at  the  same  time.  Tables  5.15  and  5.16  list  the  state 
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holding  time  distributions  for  each  of  the  states  associated  with  satellite  1  and  2, 
respectively. 


Table  5.15  State  holding  time  distributions  for  satellite  1  of  the  meteorological  constel¬ 
lation  example. 


State 

Distribution 

State 

Distribution 

1 

Gamma(5.80,0.37) 

9 

Uniform(0. 05, 1-4:0) 

2 

Gamma(4.10,6.80) 

10 

Uniform(0.13,1.10) 

3 

Gamma(1.60,2.31) 

11 

Triangular  (0.02, 0.50, 1.30) 

4 

Gamma(1.20,1.13) 

12 

Triangular(0. 13,0.40,1.20) 

5 

Erlang(2,7.41) 

13 

Triangular  (0.09, 0.30, 1.10) 

6 

Erlang(5,6.42) 

14 

Exp(1.60) 

7 

Weibull(1.20,1.22) 

15 

Exp(1.40) 

8 

Weibull(1.80,2.20) 

16 

Exp(1.80) 

Table  5.16  State  holding  time  distributions  for  satellite  2  of  the  meteorological  constel¬ 
lation  example. 


State 

Distribution 

State 

Distribution 

1 

Gamma(8.70,1.07) 

5 

Weibull(1.20,1.22) 

2 

Gamma(3.10,5.80) 

6 

Weibull(1.80,2.20) 

3 

Erlang(2,6.41) 

7 

Triangular(0. 02, 0.50, 1.30) 

4 

Erlang(5,6.42) 

8 

Exp(1.80) 
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Let  the  statistical  estimates  for  the  transition  rates  for  satellite  1  and  2  be 
given  in  the  matrices  R(1)  and  R(2\  respectively,  where 

0.064  0.013  0.019  0.019  0.013  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 

0.663-1.659  0.000  0.000  0.000  0.332  0.332  0.332  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 

0.722  0.000-1.444  0.000  0.000  0.144  0.000  0.000  0.289  0.289  0.000  0.000  0.000  0.000  0.000  0.000 

0.188  0.000  0.000-0.942  0.000  0.000  0.377  0.000  0.188  0.000  0.188  0.000  0.000  0.000  0.000  0.000 

1.482  0.000  0.000  0.000-3.705  0.000  0.000  0.371  0.000  1.482  0.371  0.000  0.000  0.000  0.000  0.000 

0.000  0.385  0.385  0.000  0.000-1.284  0.000  0.000  0.000  0.000  0.000  0.128  0.385  0.000  0.000  0.000 

0.000  0.259  0.000  0.519  0.000  0.000  -1.297  0.000  0.000  0.000  0.000  0.259  0.000  0.259  0.000  0.000 

0.000  0.742  0.000  0.000  0.742  0.000  0.000-2.474  0.000  0.000  0.000  0.000  0.742  0.247  0.000  0.000 

0.000  0.000  0.414  0.414  0.000  0.000  0.000  0.000-1.379  0.000  0.000  0.414  0.000  0.000  0.138  0.000 

0.000  0.000  0.325  0.000  0.488  0.000  0.000  0.000  0.000-1.626  0.000  0.000  0.488  0.000  0.325  0.000 

0.000  0.000  0.000  0.495  0.165  0.000  0.000  0.000  0.000  0.000-1.648  0.000  0.000  0.495  0.495  0.000 

0.000  0.000  0.000  0.000  0.000  0.347  0.520  0.000  0.347  0.000  0.000-1.734  0.000  0.000  0.000  0.520 

0.000  0.000  0.000  0.000  0.000  0.604  0.000  0.604  0.000  0.201  0.000  0.000-2.013  0.000  0.000  0.604 

0.000  0.000  0.000  0.000  0.000  0.000  0.320  0.640  0.000  0.000  0.480  0.000  0.000-1.600  0.000  0.160 

0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.280  0.420  0.420  0.000  0.000  0.000-1.400  0.280 

0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.720  0.720  0.180  0.180-1.800 

and 

-0.123  0.061  0.037  0.025  0.000  0.000  0.000  0.000 

0.748-1.871  0.000  0.000  0.561  0.561  0.000  0.000 

0.962  0.000-3.205  0.000  0.962  0.000  1.282  0.000 

0.385  0.000  0.000-1.284  0.000  0.642  0.257  0.000 
0.000  0.648  0.519  0.000-1.297  0.000  0.000  0.130 
0.000  0.990  0.000  0.742  0.000-2.474  0.000  0.742 
0.000  0.000  0.659  0.495  0.000  0.000-1.648  0.495 

0.000  0.000  0.000  0.000  0.540  0.540  0.720-1.800 

Figure  5.5  shows  the  instantaneous  availability  and  limiting  availability  for 
both  of  the  satellites  in  the  meteorological  constellation.  The  simulated  availability 
measures  were  not  included  in  this  figure  as  they  were  indistinguishable  from  the 
analytical  measures.  Table  5.17  summarizes  the  instantaneous  availability  measures 
for  each  of  the  satellites  at  t  —  10  years.  The  average  time  for  the  analytical  availabil¬ 
ity  computation  was  0.37  minutes  versus  13.03  minutes  for  the  simulation  method. 
Thus,  the  computation  increased  by  a  factor  of  over  34  for  the  simulation  results. 
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Figure  5.5  Comparison  of  instantaneous  availability  (regular)  with  the 
limiting  availability  (bold)  for  satellite  1  (regular)  and  satellite 
2  (dotted)  of  the  meteorological  constellation  example. 


Table  5.17  Sample  availability  measures  for  the  meteorological  constellation  example 


evaluated  at  t 

=  10 

years. 

Satellite 

t 

m 

Analytical  Simulated 

Difference 

1 

10 

0.885548 

0.885096 

0.000452 

2 

10 

0.660505 

0.660400 

0.000105 

The  limiting  availability  is  calculated  for  satellites  1  and  2  using  the  methods 
described  in  Chapter  4.  Table  5.18  shows  the  differences  between  the  instantaneous 
availability  and  the  limiting  availability.  Notice  that  the  limiting  availability  measure 
for  satellite  1  underestimates  the  true  availability  of  the  satellite  by  13%,  whereas 
the  limiting  availability  measure  for  satellite  2  overestimates  the  true  availability  of 
the  satellite  by  5%. 
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Table  5.18  Comparison  of  instantaneous  and  limiting  availability  measures  for  the  me¬ 
teorological  constellation  example. 


Satellite 

t 

Instantaneous 

Availability 

Limiting 

Availability 

Difference 

1 

10 

0.885548 

0.746574 

0.138974 

2 

10 

0.660505 

0.711398 

0.050893 

5-4-2  Value  Scores 

The  value  scores  of  the  meteorological  constellation  will  be  calculated  using 
the  same  set  of  attributes  that  were  used  to  score  the  satellites  of  the  previous 
constellation  examples.  The  value  functions  and  attribute  weights  for  the  three 
attributes  are  given  in  Table  5.19,  and  Table  5.20  summarizes  the  results  associated 
with  the  value  scores  for  each  of  the  satellites. 

Table  5.19  Value  functions  and  attribute  weights  for  the  meteorological  constellation 
example. 


Attribute  (i) 

y'i 

Ci 

Vi{yi) 

■uy(lO) 

1 

0.70 

1.801072 

1— exp(l. 801072 j/j) 
l—exp(l. 801072) 

0.50 

2 

0.65 

1.278652 

1 — exp(  1 .278652j/j ) 
l—exp(l. 278652) 

0.32 

3 

0.45 

-0.402692 

1— exp(— 0.402692j/j) 
l—exp(— 0.402692) 

0.18 

Table  5.20  Value  score  results  for  the  meteorological  constellation  example. 


Satellite 

Vi 

V-2 

V3 

Mv  i) 

My  2) 

My 3) 

V(10) 

1 

0.99 

0.77 

0.86 

0.978620 

0.646903 

0.883036 

0.855265 

2 

0.97 

0.77 

0.89 

0.936999 

0.646903 

0.908657 

0.839067 

5-4-3  Functional  Capability 

Based  on  the  instantaneous  availability  measurements  along  with  the  value 
scores  for  each  satellite  in  the  constellation,  the  constellation  functional  capabil¬ 
ity  measure  $(10)  is  calculated  to  be  0.655793  (0.655556  using  simulation  results). 
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Using  the  limiting  availability  measures  to  calculate  the  constellation  functional  ca¬ 
pability  measure  yields  0.617714,  which  underestimates  the  constellations  functional 
capability  by  3.8%.  An  inaccurate  constellation  functional  capability  measure  could 
lead  to  the  unnecessary  repair  or  replacement  of  satellites  in  the  constellation  in  the 
attempt  to  maintain  the  constellation’s  level  of  capability. 

5.5  Discussion 

This  chapter  illustrated  the  techniques  developed  in  Chapters  3  and  4  in  three 
different  constellations.  For  each  of  the  examples,  the  satellite  instantaneous  avail¬ 
ability  and  value  score  were  calculated  to  generate  the  constellation’s  functional  ca¬ 
pability  measure.  The  availability  measures  were  computed  using  analytical  methods 
and  Monte  Carlo  simulation,  yielding  indistinguishable  results.  This  indicates  that 
errors  associated  with  the  numerical  approximations  used  in  the  analytical  method, 
such  as  transform  approximations  and  numerical  inversions,  are  insignificant  when 
calculating  the  constellation’s  functional  capability.  Moreover,  the  analytical  results 
were  computed  with  substantially  less  effort  (i.e.,  reduced  computational  time)  than 
those  generated  via  Monte  Carlo  simulation.  A  summary  of  the  processing  times  for 
each  example  is  presented  in  Table  5.21.  The  benefits  associated  with  the  shorter 
processing  time  of  the  analytical  method  would  be  realized  for  constellations  modeled 
with  more  functions,  equating  to  a  larger  state  space. 

Table  5.21  Comparison  of  average  processing  time  (mins)  for  instantaneous  availability 
measures  for  the  constellation  examples. 


Example 

Analytical 

Simulation 

Factor  Increase 

Navstar  Constellation 

0.01 

2.62 

243.6 

Milstar  Constellation 

1.16 

10.47 

9.1 

Meteorological  Constellation 

0.37 

13.03 

34.8 

As  stated  in  Chapter  1,  the  purpose  of  the  constellation  functional  capabil¬ 
ity  measure  is  to  provide  insight  into  the  capability  of  the  constellation  in  terms  of 
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mission  effectiveness.  A  replenishment  policy  might  base  the  decision  to  perform 
a  satellite  repair  or  launch  a  new  satellite  based  on  a  minimum  threshold  of  con¬ 
stellation  functional  capability.  A  Type-I  error  associated  with  the  constellation 
capability  measure  would  imply  that  the  true  constellation  capability  is  higher  than 
what  the  measure  reports  it  to  be.  This  might  result  in  unnecessary  satellite  repairs 
or  replacements,  thus  increasing  the  overall  cost  of  maintaining  the  satellite  constel¬ 
lation.  A  Type-11  error  would  imply  that  the  true  constellation  capability  is  lower 
than  what  the  measure  reports  it  to  be.  This  might  result  in  the  failure  to  main¬ 
tain  the  constellation  above  the  minimum  threshold  of  functional  capability,  possibly 
jeopardizing  missions  which  rely  on  the  services  provided  by  the  constellation. 

In  the  first  two  examples,  the  limiting  availability  of  the  satellites  is  zero  since 
their  availability  depends  on  functions  that  cannot  be  repaired  once  they  have  failed. 
The  limiting  availability  computed  for  the  third  example  highlighted  errors  that  can 
occur  as  a  result  of  using  a  limiting  availability  in  place  of  instantaneous  availability. 
Had  the  constellation’s  functional  capability  been  assessed  at  an  earlier  time,  these 
errors  would  have  been  even  larger.  The  next  chapter  will  offer  conclusions  and 
recommendations  to  this  research  as  well  as  suggestions  for  further  research. 
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6.  Conclusions  and  Future  Research 


This  thesis  has  proposed  a  new  functional  capability  measure  for  a  satellite 
constellation.  The  methodology  is  consistent  with  the  definitions  of  performance 
measures  found  in  the  literature,  incorporating  the  constellation’s  availability  and 
effectiveness.  Chapter  1  introduced  constellation  capability  measures  and  discussed 
their  alignment  with  capability-based  assessment  of  U.S.  military  systems.  This 
measure  is  a  necessary  component  for  capability-based  replenishment  policies,  pos¬ 
sibly  giving  space  system  planners  information  to  justify  program  funds.  Chapter  2 
showed  how  this  research  extends  the  literature;  a  constellation  specific  measure  is 
developed  that  incorporates  a  stochastic  model  of  satellite  availability  with  multi¬ 
attribute  value  theory. 

Chapters  3  and  4  respectively  presented  mathematical  models  for  the  explicit 
computation  of  the  functional  capability  measure.  Each  satellite’s  instantaneous 
availability  and  value  score  are  calculated  at  a  common  constellation  age  yielding 
a  functional  capability  score.  The  constellation  measure  is  formed  by  averaging 
these  capability  scores.  By  employing  a  few  reasonable  assumptions,  the  satellite 
degradation  is  modeled  as  a  stochastic  process  to  acquire  its  instantaneous  avail¬ 
ability.  Finally,  Chapter  5  illustrated  the  application  of  the  methodology  in  three 
constellation  examples.  The  functional  capability  measures  were  computed  for  a 
constellation  via  individual  satellite  availability  and  value  scores.  Results  of  the  in¬ 
stantaneous  availability  computation  were  compared  with  those  obtained  by  Monte 
Carlo  simulation  and  were  shown  to  be  consistent. 

The  methodology  proposed  in  this  thesis  results  in  a  mathematically  tractable 
solution,  allowing  the  measure  to  be  directly  computed.  The  speed  of  computa¬ 
tion  for  the  instantaneous  availability  measures  was  considerably  faster  than  that 
required  for  the  Monte  Carlo  simulation  results.  The  errors  associated  with  the 
numerical  inversion  and  transform  approximation  methods  of  the  analytical  solu- 
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tion  were  negligible,  yielding  an  indistinguishable  constellation  capability  measure. 
Incorporating  the  instantaneous  availability  model  to  a  satellite  with  M  functions, 
does  require  computation  of  matrices  of  dimension  2M .  Satellites  with  many  func¬ 
tions  resulting  in  large  state  spaces  could  possibly  lead  to  computational  complexity 
issues. 

The  imposed  assumptions  for  the  stochastic  degradation  process  of  the  satellite 
are  not  entirely  unrealistic.  The  model  does  not  require  satellites  to  be  identical. 
The  meteorological  constellation  example  presented  in  Chapter  5  demonstrated  this 
flexibility.  Also,  in  the  case  of  a  satellite  with  exponential  function  lifetimes,  the 
methodology  of  Chapter  3  is  exact.  The  methodology  of  Chapter  4  to  a  satellite 
with  non-exponential  function  lifetimes  assumes  the  sojourn  time  distributions  are 
independent  of  the  next  transition  state.  This  may  be  an  unrealistic  assumption 
given  the  manor  in  which  the  satellite’s  state  space  is  constructed.  However,  a  model 
requiring  all  of  the  conditional  sojourn  time  distributions  between  all  possible  states 
would  be  infeasable,  as  acquiring  data  to  parametrically  estimate  these  distributions 
is  nearly  impossible. 

The  independence  assumption  for  functions  and  satellites  may  also  be  unre¬ 
alistic.  The  main  subsystems  of  a  satellite  are  designed  with  as  much  redundancy 
as  possible  to  prevent  catastrophic  satellite  failure  in  the  event  that  one  function 
fails.  In  reality,  the  satellite  payload  depends  on  the  satellite  bus,  which  provides  a 
platform  for  the  mission  related  functions.  In  terms  of  satellite  independence,  some 
constellations  have  designated  relay  satellites,  which  act  as  the  communication  link 
from  ground  stations  to  the  constellation.  Other  satellites  within  the  constellation 
are  dependent  on  this  relay  satellite  to  perform  their  functions.  Therefore,  this  model 
assumption  may  be  violated  for  certain  constellations. 

The  methodology  is  not  specific  to  one  type  of  constellation  mission.  The 
satellite  value  scores  are  based  on  attributes  constructed  to  directly  reflect  mission 
effectiveness  of  the  constellation.  However,  the  value  scores  may  be  entirely  sub- 
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jective.  These  scores  depend  on  accurate  information  from  a  knowledgable  subject 
matter  expert  to  construct  and  score  each  of  the  satellites  in  the  constellation.  The 
value  score  is  ultimately  the  opinion  of  a  SME;  constellation  measures  may  not  be 
consistent  when  judged  with  different  SMEs.  Moreover,  there  are  many  methods  for 
generating  a  multi-attribute  value;  this  thesis  presented  only  one  of  them.  Certain 
methods  may  be  more  advantageous  to  specific  satellite  constellations,  and  also  to 
particular  subject  matter  experts. 

Our  proposed  model  for  constellation  functional  capability  can  benefit  from 
improvements  which  generalize  and  facilitate  its  application  so  we  now  discuss  a  few 
suggestions  for  model  improvements.  First,  this  model  only  encompasses  the  space 
segment  of  an  entire  system  that  accomplishes  space-based  missions.  The  model  can 
be  expanded  to  include  the  ground  segment  of  the  system  as  well.  Second,  the  idea 
of  on-orbit  satellite  spares  could  be  incorporated  into  the  model.  A  possible  way  to 
accomplish  this  might  involve  constructing  the  state  space  in  a  manner  which  allows 
multiple  failed  functions  to  be  brought  back  to  an  available  state  simultaneously, 
simulating  the  entire  replacement  of  a  satellite. 

The  method  for  constructing  the  satellite  value  scores  is  another  area  that  can 
be  improved.  In  particular,  the  value  score  may  be  incorporated  as  a  reward  within 
a  Markov  reward  model  that  can  be  analyzed  using  a  number  of  existing  techniques. 
Finally,  a  different  approach  to  modeling  satellites  with  non-exponential  function 
lifetimes  is  to  use  phase-type  distributions,  allowing  non-exponential  distributions 
to  be  transformed  into  exponential  distributions. 
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Appendix  A.  Computer  Code 

A.l  Analytical  Availability 

function  [t,A_SMP,A_stat, duration]  =  Availability_Examplel (t) 

7.  Author:  Capt  Cole  Gulyas 

7.  Date:  19  Jan  05 

7.  Input:  t=time 

7,  Output:  Measures  for  Example  1:  t=time,  A_SMP=instantaneous  availability, 

7,  A_stat=stationary  availability,  and  duration=processing  time. 

7.  Subfunctions:  prob_matrix(Q , t)  by  Capt  Cole  Gulyas,  Dr.  Jeffrey  Kharoufeh 
7,  invt_lap(Q ,t ,row, col)  by  Dr.  Jeffrey  Kharoufeh 

7,  el35a(Q ,x,y ,row, col)  by  Dr.  Jeffrey  Kharoufeh,  Capt  Cole  Gulyas 

7,  rate_matrix(m,f  ,r)  by  Capt  Cole  Gulyas 

tic  7.  Start  timer 

71-4  Block  A 
m=3 ; 
n=2~m ; 

a=[l  0  0  0  0  0  0  0]  ; 
s=[l  0  0  0  0  0  0  0]  ; 
f =  [ . 0291  .0503  .0230]; 
r= [0  3.3300  0] ; 

7,  7.5-16  Block  A 
7.  m=3; 

7.  n=2~m ; 

7.  a=[l  0  0  0  0  0  0  0]  ; 

7.  s=[l  0  0  0  0  0  0  0]  ; 

7.  f  =  [  .0121  .0421  .0189]  ; 

7.  r=  [0  4.8216  0]  ; 

7.  7.17-20  Block  R 
7#  m=3; 

7.  n=2~m; 

7.  a=  [ .  98  0  . 02  0  0  0  0  0]  ; 

7.  s=[l  0  0  0  0  0  0  0]  ; 

7.  f  =  [  .0121  .0303  .0143]  ; 

7.  r=  [0  4.9216  0]  ; 

7.  7.21-24  Block  R 
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7,  m=3; 

7.  n=2''m; 

7.  a=[l  0  0  0  0  0  0  0]  ; 

'/,  s=[l  0  0  0  0  0  0  0]; 

7.  f  =  [  .0103  .0282  .0118]  ; 

7.  r=  [0  4.9216  0]  ; 

Q=rate_matrix(m, f ,r) ; 

7,  Mean  holding  times  for  each  state 
mu=-l . /diag(Q) ; 

7.  P*(s) 

P_Laplace_SMP=prob_matrix(Q ,t) ; 

7,  Pmf 

pmf _SMP=a*P_Laplace_SMP ; 

7,  Availability 
A_SMP=pmf _SMP*s ’ ; 

duration=toc ;  7«  Stop  timer 

7,  Calculate  stationary  distribution  and  availability 
7.  Pi 

ze=zeros(l ,n) ; 
ze(end)=l ; 

Pe=  []  ; 
for  i=l:n 

for  j=l:n 

if  i~=j  Pe(i, j)=-Q(i, j)/Q(i,i) ;  end 

end 

end 

Pe=Pe-eye (n) ; 

Pe( : , end)=ones (n, 1) ; 
pie=ze*inv(Pe) ; 

7,  Pmf  and  stationary  availability 

p=n ; 

for  j=l:n 

p=[p  (pie(j)*mu(j))/(pie*(mu))] ; 

end 

A_stat=p*s  * ; 
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7, - Sub  Function - 

function  eq  =  el35a(Q ,x,y ,row, col) 

7,  Author:  Capt  Gulyas,  Dr.  Jeffrey  Kharoufeh 


7.  Last  Revision:  4  Oct  04 

7.  Note:  function  handles  Q,  row,  and  col  were  added  to  function  el35a 

s=x+y*i;  7*  Complex  transform  variable 

P_star  =  []  ;  7*  Declare  matrix 

I=eye(size(Q)) ;  %  Identity  matrix 

7.  All  states  are  distributed  exponential: 

LSTG=  []  ; 

LSTD=  []  ; 

for  a=l : length(Q) 

for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*(-Q(a,a)/(s-Q(a,a))) ;  end 
LSTD(a,a)=s/(s-Q(a,a)) ; 

end 

end 

P_star  =  (1/s)  *  (inv(I-LSTG)  )  *LSTD;  7.  Build  P*(s) 

eq  =  real (P_star (row, col) ) ;  7»  Return  first  moment 

7. - End  Sub  Function - 

function  [t ,A_SMP,A_stat , duration]  =  Availability_Example2 (t) 

7.  Author:  Capt  Cole  Gulyas 

7.  Date:  5  Jan  05 

7.  Input:  t=time 

7.  Output:  Measures  for  Example  2:  t=time,  A_SMP=instantaneous  availability, 

7»  A_stat=stat ionary  availability,  and  duration=processing  time. 

7.  Subfunctions:  prob_matrix(Q , t)  by  Capt  Cole  Gulyas,  Dr.  Jeffrey  Kharoufeh 
7.  invt_lap(Q ,t ,row, col)  by  Dr.  Jeffrey  Kharoufeh 

7.  el35a(Q,x,y,row,col)  by  Dr.  Jeffrey  Kharoufeh,  Capt  Cole  Gulyas 

tic  7.  Start  timer 

t=10 ; 
m=4 ; 
n=2~m; 

a=  [1  00000000000000  0]; 
s=  [1  01110000000000  0]; 

K=  [3.8  .37  0  7o  Gamma 
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2.1 

.50 

0 

°/0  Gamma 

2.2 

.41 

0 

°/0  Gamma 

3.2 

.42 

0 

7,  Gamma 

1.2 

.22 

0 

7.  Weibull 

1.8 

.42 

0 

7,  Weibull 

1.9 

.31 

0 

7,  Weibull 

1.2 

.33 

0 

7.  Weibull 

1.8 

2.1 

2.3 

7,  Triangular 

1.63 

.  2.2 

3.2 

7,  Triangular 

1.1 

1.83 

2.4 

7,  Triangular 

1.7 

2.5 

4.51 

7,  Triangular 

.05 

2.3 

0 

7o  Uniform 

.  13 

2.1 

0 

7,  Uniform 

.03 

2.4 

0 

7,  Uniform 

.51 

.93 

0]; 

7o  Uniform 

'/,  Mean  holding  times  for  each  state 
mu=  []  ;  for  i=l  :4 

mu=[mu;  K(i ,  1)  /K(i  ,  2)]  ; 
end  for  i=5:8 

mu= [mu;  (1/K(i , 2) ) *gamma( (K(i ,1)+1)/K(i,l))] ; 
end  for  i=9:12 


"/.Gamma 

•/.Weibull 


mu=[mu;  (K(i ,  1)+K(i ,2)+K(i  , 3) ) /3]  ;  "/.Triangular 

end  for  i=13:16 


mu=[mu;  (K(i ,  1)+K(i  ,2)  ) /2]  ;  "/.Uniform 

end 


"/,  Rate  matrix  built  with  only  function  1  repaired  if  function  2  is 
"/,  available 

R=[-l/mu(l)  . 2/mu(l)  ,3/mu(l)  .3/mu(l)  ,2/mu(l)  00000000000 
. 2/mu (2)  -l/mu(2)  000  .4/mu(2)  ,2/mu(2)  .2/mu(2)  00000000 
0  0  -l/mu(3)  0  0  . 2/mu(3)  0  0  . 4/mu (3)  .4/mu(3)  000000 
000  -l/mu(4)  0  0  .5/mu (4)  0  . 3/mu (4)  0  . 2/mu (4)  00000 
0000  -l/mu(5)  0  0  ,2/mu(5)  0  .7/mu(5)  . 1/mu (5)  00000 
0  0  . 4/mu(6)  0  0  -l/mu(6)  00000  .3/mu(6)  .3/mu(6)  000 
000  . 5/mu (7)  0  0  -l/mu(7)  0000  ,2/mu(7)  0  ,3/mu(7)  0  0 
0000  .4/mu (8)  0  0  -l/mu(8)  0000  .5/mu(8)  ,l/mu(8)  0  0 
00000000  -l/mu(9)  0  0  ,7/mu(9)  0  0  .3/mu(9)  0 
000000000  -l/mu(10)  0  0  ,8/mu(10)  0  .2/mu(10)  0 
0000000000  -l/mu(ll)  0  0  ,6/mu(ll)  .4/mu(ll)  0 
00000000  . 2/mu(12)  0  0  -l/mu(12)  000  .8/mu(12) 
000000000  .3/mu (13)  0  0  -l/mu(13)  0  0  ,7/mu(13) 
0000000000  .3/mu (14)  0  0  -l/mu(14)  0  .7/mu(14) 
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00000000000000  -l/mu(15)  l/mu(15) 
000000000000000  eps] ; 

7.  P*(s) 

P_Laplace_SMP=prob_matrix(R,t) ; 

7,  Pmf 

pmf _SMP=a*P_Laplace_SMP ; 

7*  Availability 
A_SMP=pmf _SMP*s ’ ; 

duration=toc  l  Stop  timer 

7,  Calculate  stationary  distribution  and  availability 
7.  Pi 

ze=zeros(l ,n) ; 
ze(end)=l ; 

Pe=  []  ; 
for  i=l:n 

for  j=l:n 

if  i~=j  Pe(i, j)=-R(i, j)/R(i,i) ;  end 

end 

end 

Pe=Pe-eye (n) ; 

Pe( : , end)=ones (n, 1) ; 
pie=ze*inv(Pe) ; 

7,  Pmf  and  stationary  availability 
p=[]  ;  for  j=l  :n 

p=[p  (pie(j)*mu(j))/ (pie* (mu))] ; 
end  A_stat=p*s’; 

7, - Sub  Function - 

function  eq  =  el35a(Q ,x,y ,row, col) 

7,  Author:  Capt  Gulyas,  Dr.  Jeffrey  Kharoufeh 

7.  Last  Revision:  4  Oct  04 

7,  Note:  function  handles  Q,  row,  and  col  were  added  to  function  el35a 

s=x+y*i;  %  Complex  transform  variable 

P_star  =  [] ;  %  Declare  matrix 

I=eye(size(Q)) ;  %  Identity  matrix 
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K=  [3.8 

.37 

0 

7.  Gamma 

2.1 

.50 

0 

7.  Gamma 

2.2 

.41 

0 

7,  Gamma 

3.2 

.42 

0 

7.  Gamma 

1.2 

.22 

0 

7.  Weibull 

1.8 

.42 

0 

7.  Weibull 

1.9 

.31 

0 

7.  Weibull 

1.2 

.33 

0 

7.  Weibull 

1.8 

2.1 

2.3 

7,  Triangular 

1.61 

.  2.2 

3.2 

7,  Triangular 

1.1 

1.83 

2.4 

7.  Triangular 

1.7 

2.5 

4.51 

7,  Triangular 

.05 

2.3 

0 

7,  Uniform 

.  13 

2.1 

0 

7,  Uniform 

.03 

2.4 

0 

7,  Uniform 

.51 

.93 

0]; 

7.  Uniform 

LSTG=  []  ;  LSTD=  []  ; 

7,  States  1-4:  Gamma 
for  a=l:4 

F_lst=(K(a,2)/(s+K(a,2)))~K(a,l)  ; 

for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst ; 

end 

end 

7*  States  5-8:  Weibull 
for  a=5:8 

sumterm=0; 

N=1000 ; 
for  k=l : N 

sumterm=sumterm  +  exp(-s* (1/K(a, 2) ) * (-log(l- (k/(N+l) ) ) ) ~ (1/K(a, 1) ) ) ; 

end 

F_lst=(l/N) *sumterm; 
for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst; 

end 

end 

7,  States  9-12:  Triangular 


A-6 


for  a=9:12 


terml=(exp(-s*K(a, 1) )-exp(-s*K(a, 2) ) ) /( (K(a,3) -K(a, 1) ) * (K(a,2) -K(a, 1) ) ) ; 
term2=(exp(-s*K(a,3) )-exp(-s*K(a, 2) ) ) /( (K(a,3) -K(a, 1) ) * (K(a,3) -K(a,2) ) ) ; 
F_lst=(2/(s~2))*(terml+term2)  ; 
for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst;  end 
LSTD(a,a)=l-F_lst ; 

end 

end 

7,  States  13-16:  Uniform 
for  a=13 : 16 

F_lst=(exp(-s*K(a, 1) )-exp(-s*K(a, 2) ) ) /(s* (K(a, 2)-K(a, 1) ) ) ; 
for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst; 

end 

end 

P_star  =  (1/s) * (inv(I-LSTG)  )*LSTD;  7.  Build  P*(s) 

eq  =  real (P_star (row, col) )  ;  7*  Return  first  moment 

7. - End  Sub  Function - 


function  [t , A_SMP ,A_stat , duration]  =  Availability_Example3_l (t) 

7,  Author:  Capt  Cole  Gulyas 

7,  Date:  22  Jan  05 

7.  Input:  t=time 

7,  Output:  Measures  for  Example  3,  satellite  1:  t=time,  A_SMP=instantaneous 

7,  availability,  A_stat=stat ionary  availability,  and  duration=pro cessing  time. 

7,  Subfunctions:  prob_matrix(Q , t)  by  Capt  Cole  Gulyas,  Dr.  Jeffrey  Kharoufeh 
7,  invt_lap(Q,t , row, col)  by  Dr.  Jeffrey  Kharoufeh 

7.  el35a(Q ,x,y ,row, col)  by  Dr.  Jeffrey  Kharoufeh,  Capt  Cole  Gulyas 


tic 

7,  Start 

timer 

m=4 ; 

n=2~m; 

o 

o 

II 

0  0 

0  0  0 

0  0  0 

0 

0  0  0  0]; 

o 

o 

II 

CQ 

0  0 

0  0  0 

0  0  0 

0 

0  0  0  0]; 

II 

on 

00 

.37  0 

7o  Gamma 

4.1 

6.80 

0 

7,  Gamma 

1.6 

2.31 

0 

7,  Gamma 

1.2 

1.13 

0 

7.  Gamma 
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2  7.41  0 

7,  Erlang 

5  6 . 42  0 

7,  Erlang 

1.2  1.22  0 

7.  Weibull 

1.8  2.2  0 

7.  Weibull 

.05  1.4  0 

7,  Uniform 

.13  1.1  0 

7,  Uniform 

.02  .5  1.3 

7,  Triangular 

.13  .4  1.2 

7,  Triangular 

.09  .3  1.1 

7,  Triangular 

11.6  0 

7«  Exp 

11.4  0 

7.  Exp 

1  1.8  0  ] ; 

7o  Exp 

ian  holding  times 

for  each  state 

mu=  []  ; 

for  i=l:6 

mu= [mu;  K(i , 1) /K(i ,2)] ; 

end 

for  i=7:8 

mu= [mu;  (1/K(i ,2) ) *gamma( (K(i , 1)+1)/K(i , 1) )] ; 

end 

for  i=9:10 

mu= [mu;  (K(i , 1) +K(i , 2) ) / 2]  ; 

end 

for  i=ll:13 

mu= [mu;  (K(i , 1) +K(i , 2)+K(i , 3) ) / 3]  ; 

end 

for  i=14:16 

mu= [mu;  K(i , 1) /K(i ,2)] ; 

end 


°/oGamma  &  Erlang 


•/.Weibull 


°/0Unif  orm 


70Tri  angular 


°/0Exp 


7.  Rate  matrix  built  with  all  types  of  repairs  possible 
R=[-l/mu(l)  . 2/mu(l)  .3/mu(l)  .3/mu(l)  .2/mu(l)  00000000000 
. 4/mu (2)  -l/mu(2)  000  .2/mu(2)  .2/mu(2)  .2/mu(2)  00000000 
. 5/mu (3)  0  -l/mu(3)  0  0  .l/mu(3)  0  0  .2/mu(3)  .2/mu(3)  000000 
. 2/mu (4)  0  0  -l/mu(4)  0  0  .4/mu(4)  0  .2/mu(4)  0  .2/mu(4)  00000 
.4/mu(5)  000  -l/mu(5)  0  0  .l/mu(5)  0  .4/mu(5)  .l/mu(5)  00000 
0  . 3/mu(6)  . 3/mu(6)  0  0  -l/mu(6)  00000  .l/mu(6)  .3/mu(6)  000 
0  . 2/mu(7)  0  . 4/mu (7)  0  0  -l/mu(7)  0000  .2/mu(7)  0  .2/mu(7)  0  0 
0  . 3/mu(8)  0  0  . 3/mu(8)  0  0  -l/mu(8)  0000  .3/mu(8)  .l/mu(8)  0  0 
0  0  . 3/mu(9)  . 3/mu(9)  0000  -l/mu(9)  0  0  .3/mu(9)  0  0  .l/mu(9)  0 
0  0  . 2/mu(10)  0  . 3/mu(10)  0000  -l/mu(10)  0  0  .3/mu(10)  0  .2/mu(10)  0 
000  . 3/mu(ll)  . l/mu(ll)  00000  -l/mu(ll)  0  0  .3/mu(ll)  .3/mu(ll)  0 
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00000  . 2/mu(12)  .3/mu(12)  0  .2/mu(12)  0  0  -l/mu(12)  000 
00000  .3/mu (13)  0  .3/mu(13)  0  .l/mu(13)  0  0  -l/mu(13)  0  0 
000000  .2/mu (14)  .4/mu(14)  0  0  .3/mu (14)  0  0  -l/mu(14)  0 
00000000  . 2/mu(15)  .3/mu(15)  .3/mu(15)  000  -l/mu(15) 
00000000000  . 4/mu(16)  .4/mu(16)  .l/mu(16)  .l/mu(16) 


.3/mu (12) 
.3/mu (13) 

. l/mu(14) 

. 2/mu(15) 
-l/mu(16)] ; 


1  P*(s) 

P_Laplace_SMP=prob_matrix(R,t) ; 


7,  Pmf 

pmf _SMP=a*P_Laplace_SMP ; 


7#  Availability 
A_SMP=pmf _SMP*s ’ ; 

duration=toc ;  7*  Stop  timer 

'/.  Calculate  stationary  distribution  and  availability 

'/.  Pi 

ze=zeros(l ,n) ; 
ze(end)=l ; 

Pe=  []  ; 
for  i=l:n 

for  j=l:n 

if  i~=j  Pe(i, j)=-R(i, j)/R(i,i) ;  end 

end 

end 

Pe=Pe-eye (n) ; 

Pe( : , end)=ones (n, 1) ; 
pie=ze*inv(Pe) ; 

'/.  Pmf  and  stationary  availability 

p=n ; 

for  j=l:n 

p=[p  (pie(j)*mu(j))/(pie*(mu))] ; 

end 

A_stat=p*s ’ ; 


°/0 - Sub  Function - 

function  eq  =  el35a(Q ,x,y ,row, col) 

7#  Author:  Capt  Gulyas,  Dr.  Jeffrey  Kharoufeh 


7,  Last  Revision:  4  Oct  04 

%  Note:  function  handles  Q,  row,  and  col  were  added  to  function  el35a 
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s=x+y*i ; 

P_star  =  []  ; 
I=eye(size(Q)) ; 


[5.8  .37  0 

7o  Gamma 

4.1  6.80  0 

7.  Gamma 

1.6  2.31  0 

7.  Gamma 

1.2  1.13  0 

7.  Gamma 

2  7.41  0 

7.  Erlang 

5  6 . 42  0 

7.  Erlang 

1.2  1.22  0 

7.  Weibull 

1.8  2.2  0 

7.  Weibull 

.05  1.4  0 

7,  Uniform 

.13  1.1  0 

7.  Uniform 

.02  .5  1.3 

7.  Triangular 

.13  .4  1.2 

7.  Triangular 

.09  .3  1.1 

7.  Triangular 

11.6  0 

7»  Exp 

11.4  0 

7«  Exp 

1  1.8  0  ] ; 

7*  Exp 

LSTG=  []  ; 

LSTD=  []  ; 

7,  States  1-6:  Gamma  and  Erlang 
for  a=l:6 

F_lst=(K(a,2)/(s+K(a,2)))~K(a,l)  ; 

for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst; 

end 

end 

7.  States  7-8:  Weibull 
for  a=7:8 

sumterm=0 ; 

N=1000; 
for  k=l:N 

sumterm=sumterm  +  exp(-s* (1/K(a, 2) ) * (-log(l- (k/(N+l) ) ) ) ~ (1/K(a, 1) ) ) ; 

end 

F_lst=(l/N) *sumterm; 
for  b=l : length(Q) 


/  Complex  transform  variable 
7,  Declare  matrix 
7,  Identity  matrix 
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if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst; 

end 

end 

7,  States  9-10:  Uniform 
for  a=9:10 

F_lst=(exp(-s*K(a, 1) )-exp(-s*K(a, 2) ) ) /(s* (K(a, 2)-K(a, 1) ) ) ; 
for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst; 

end 

end 

7*  States  11-13:  Triangular 
for  a=ll:13 

terml=(exp(-s*K(a, 1) )-exp(-s*K(a, 2) ) ) /( (K(a,3) -K(a, 1) ) * (K(a,2) -K(a, 1) ) ) ; 
term2=(exp(-s*K(a,3) )-exp(-s*K(a, 2) ) ) /( (K(a,3) -K(a, 1) ) * (K(a,3) -K(a,2) ) ) ; 
F_lst=(2/(s'‘2))*(terml+term2)  ; 
for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst;  end 
LSTD(a,a)=l-F_lst ; 

end 

end 

7,  States  14-16:  Exponential 
for  a=14:16 

F_lst=(K(a,2)/  (s+K(a,2)  )  )  "'K(a,  1)  ; 
for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst; 

end 

end 

P_star  =  (1/s)  *  (inv(I-LSTG)  )  *LSTD;  7.  Build  P*(s) 

eq  =  real (P_star (row, col) )  ;  7*  Return  first  moment 

7, - End  Sub  Function - 


function 
7,  Author: 
7,  Date: 

7.  Input : 


[t ,  A_SMP ,  A_stat ,  duration] 
Capt  Cole  Gulyas 
22  Jan  05 
t=time 


Availability_Example3_2 (t) 


A-ll 


7,  Output:  Measures  for  Example  3,  satellite  2:  t=time,  A_SMP=instantaneous 

7*  availability,  A_stat=stationary  availability,  and  duration=pro cessing  time. 

7,  Subfunctions:  prob_matrix(Q , t)  by  Capt  Cole  Gulyas,  Dr.  Jeffrey  Kharoufeh 
7,  invt_lap(Q ,t ,row, col)  by  Dr.  Jeffrey  Kharoufeh 

7,  el35a(Q,x,y,row,col)  by  Dr.  Jeffrey  Kharoufeh,  Capt  Cole  Gulyas 

tic  7.  Start  timer 

m=3 ; 
n=2~m; 

a=[l  0  0  0  0  0  0  0]  ; 
s=[l  0  0  0  0  0  0  0]  ; 


[8.7  1.07  0 

7o  Gamma 

3.1  5.80  0 

7.  Gamma 

2  6.41  0 

7.  Erlang 

5  6 . 42  0 

7.  Erlang 

1.2  1.22  0 

7.  Weibull 

1.8  2.2  0 

7.  Weibull 

.02  .5  1.3 

7.  Triangular 

1  1.8  0  ] ; 

7*  Exp 

7,  Mean  holding  times  for  each  state 
mu=  []  ; 
for  i=l:4 

mu=  [mu;  K(i,  1)/K(i,2)]  ;  7oGamma  &  Erlang 

end 

for  i=5:6 

mu=  [mu;  (1/K(i ,2) ) *gamma( (K(i, 1)+1)/K(i , 1) )] ;  7oWeibull 

end 

for  i=7 

mu=  [mu;  (K(i,l)+K(i,2)+K(i,3))/3]  ;  7»Triangular 

end 

for  i=8 

mu=  [mu;  K(i,l)/K(i,2)]  ;  7«Exp 

end 

7.  Rate  matrix  built  with  all  types  of  repairs  possible 
7,  R  matrix  (observed  rates  of  state  transitions) 

R=[  -l/mu(l)  . 5/mu(l)  .3/mu(l)  .2/mu(l)  0000 
.4/mu(2)  -l/mu(2)  0  0  .3/mu(2)  .3/mu(2)  0  0 

. 3/mu(3)  0  -l/mu(3)  0  .3/mu(3)  0  .4/mu(3)  0 

.3/mu(4)  0  0  -l/mu(4)  0  .5/mu(4)  .2/mu(4)  0 

0  .5/mu(5)  .4/mu(5)  0  -l/mu(5)  0  0  .l/mu(5) 
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0  .4/mu(6)  0  .3/mu(6)  0  -l/mu(6)  0  .3/mu(6) 

0  0  .4/mu (7)  . 3/mu (7)  0  0  -l/mu(7)  .3/mu(7) 
0000  .3/mu(8)  .3/mu(8)  .4/mu(8)  -l/mu(8)] ; 

7.  P*(s) 

P_Laplace_SMP=prob_matr ix (R , t ) 

7.  Pmf 

pmf _SMP=a*P_Laplace_SMP 

°/«  Availability 
A_SMP=pmf _SMP*s ’ 

duration=toc ;  °/e  Stop  timer 

°/«  Calculate  stationary  distribution  and  availability 

7.  Pi 

ze=zeros(l ,n) ; 
ze(end)=l ; 

Pe=  []  ; 
for  i=l:n 

for  j=l:n 

if  i~=j  Pe(i, j)=-R(i, j)/R(i,i) ;  end 

end 

end 

Pe=Pe-eye (n) ; 

Pe( : , end)=ones (n, 1) ; 
pie=ze*inv(Pe) ; 

°/0  Pmf  and  stationary  availability 

p=[] ; 
for  j=l:n 

p=[p  (pie(j)*mu(j))/ (pie* (mu))] ; 

end 

A_stat=p*s * ; 


7, - Sub  Function - 

function  eq  =  el35a(Q ,x,y ,row, col) 

7.  Author:  Capt  Gulyas,  Dr.  Jeffrey  Kharoufeh 


7.  Last  Revision:  4  Oct  04 

7.  Note:  function  handles  Q,  row,  and  col  were  added  to  function  e!35a 


s=x+y*i;  %  Complex  transform  variable 
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P_star  =  []  ; 
I=eye(size(Q)) ; 


[8.7  1.07  0 

7.  Gamma 

3.1  5.80  0 

7.  Gamma 

2  6.41  0 

7,  Erlang 

5  6 . 42  0 

7,  Erlang 

1.2  1.22  0 

7.  Weibull 

1.8  2.2  0 

7.  Weibull 

.02  .5  1.3 

l  Triangular 

1  1.8  0  ] ; 

7*  Exp 

°/e  Declare  matrix 
°/0  Identity  matrix 


LSTG=  []  ; 
LSTD=  []  ; 


7,  States  1-4:  Gamma  and  Erlang 
for  a=l:4 

F_lst=(K(a,2)/(s+K(a,2)))~K(a,l)  ; 

for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst ; 

end 

end 


7,  States  5-6:  Weibull 
for  a=5:6 

sumterm=0; 

N=1000; 
for  k=l : N 

sumterm=sumterm  +  exp(-s* (1/K(a, 2) ) * (-log(l- (k/(N+l) ) ) ) ~ (1/K(a, 1) ) ) ; 

end 

F_lst=(l/N) *sumterm; 
for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst ; 

end 

end 


7,  State  7 :  Triangular 
for  a=7 

terml=(exp(-s*K(a, 1) )-exp(-s*K(a, 2) ) ) /( (K(a,3)-K(a, 1) ) * (K(a,2) -K(a, 1) ) ) ; 
term2=(exp(-s*K(a,3) )-exp(-s*K(a, 2) ) ) /( (K(a,3)-K(a, 1) ) * (K(a,3)-K(a,2) ) ) ; 
F_lst=(2/(s~2))*(terml+term2) ; 
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for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst; 

end 

end 

7,  State  8:  Exponential 
for  a=8 

F_lst=(K(a,2)/  (s+K(a,2)  )) ~K(a, 1)  ; 
for  b=l : length(Q) 

if  a~=b  LSTG(a,b)=(Q(a,b)/(-Q(a,a)))*F_lst ;  end 
LSTD(a,a)=l-F_lst; 

end 

end 

P_star  =  (1/s) * (inv(I-LSTG) )*LSTD;  7.  Build  P*(s) 

eq  =  real (P_star (row, col) )  ;  7«  Return  first  moment 

7, - End  Sub  Function - 


function  fl  =  invt_lap(Q ,t ,row, col) 

7,  The  purpose  of  this  MATLAB  program  is  to  approximate  the  inverse  transform  of  a  one- 
7.  dimensional  Laplace  transform  in  order  to  find  the  moments  of  the  probability 
7.  distribution,  G(t) .  The  program  is  based  on  the  algorithm  of  Abate  and  Whitt  (1995) . 


7. 


7,  Author : 
7,  Date : 
7.  Last  Revised: 
7,  References: 
7. 


Jeffrey  P.  Kharoufeh,  Ph.D.  Candidate,  IE  &  OR,  Penn  State  University 
January  23,  2001 
February  5,  2001 

Abate,  J.  and  W.  Whitt  (1995).  Numerical  Inversion  of  the  Laplace 
Transform  of  Probability  Distribution.  ORSA  Journal  on  Computing,  7, 


7. 


36-43. 


7,  Note:  function  handles  Q,  row,  and  col  were  added  to  function  invt_lap 
7«Initialize  variables,  set  parameters 

rho=0.8;  qx=[0.8];  tx=[0];  m-11;  c=[];  ga=8;  A=ga*log(10) ;  mm=2'‘m; 

for  k=0:m 

d=nchoosek(m,k) ; 
c=  [c  d]  ; 

end 

for  t  =  t; 
tx  =  t ; 
ntr=15 ; 
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u=exp(A/2)/t ; 
x=A/ (2*t) ; 
h=pi/t ; 

su=zeros (m+2) ; 
sm=e 135a (Q ,x,0 , row, col) /2 ; 
for  k=l:ntr 
y=k*h; 

sm=sm+((-l) ~k)*el35a(Q,x,y,row,col) ; 

end 

su(l)=sm; 
for  k=l : 12 
n=ntr+k; 
y=n*h; 

su(k+l)=su(k)+( (-1) ~n) *e 135a (Q ,x,y, row, col) ; 

end 

avl=0;  av2=0; 
for  k=l : 12 

avl=avl+c(k)*su(k) ; 
av2=av2+c(k)*su(k+l) ; 

end 

fl  =  u*avl/min;  f 2=u*av2/mm;  qx=[qx  f 2]  ; 

end 


function  P  =  prob_matrix(Q , t) 

7,  Authors:  Capt  Cole  Gulyas,  Dr.  Jeffrey  Kharoufeh 

7,  Date:  27  Oct  04 

7o  Input:  Q=rate  transition  matrix,  t=Time,  and  eventually  state 

7,  distributions,  etc. 

7.  Output:  Approximation  of  P(t),  the  probability  transition  matrix,  using 

7,  functions  INVT_LAP  and  el35  created  by  Dr.  Kharoufeh  as  subfunctions. 

7.  References:  Kharoufeh,  Jeffrey.  INVT_LAP.  MatLab  code.  Feb  5,  2001. 

7.  Kharoufeh,  Jeffrey.  el35.m.  MatLab  code.  Sept.  22,  2004. 

7,  See  references  for  subfunctions. 


format  long 


7.  Declare  matrix 
P=[]  ; 


7oBuild  P  matrix  elements  (required  since  INVT_LAP  subfunction  is  univariate) 
for  i  =  l:length(Q) 
for  j=l : length(Q) 

disp([’(J  num2str(i)  3  , 3  num2str(j)  ’)’]) 
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P(i, j)=invt_lap(Q ,t ,i, j) ; 


end 

end 


function  Q 
7,  Author: 

7,  Date: 

7.  Input : 

7. 

7.  Output : 


rate_matrix(m,f ,r) 

Capt  Cole  Gulyas 
25  Sep  04 

m=number  of  functions  in  a  satellite,  f=lxm  vector  of 
function  failure  rates,  r=lxm  vector  of  repair  rates. 
Rate  matrix  Q. 


7,  Initialize  the  storage  array  for  S. 
S=[]  ; 


7,  Determine  dimension  of  S. 
dimS=2''m; 


7.  Create  a  matrix  were  the  rows  will  be  permuted  to  represent  the 
7,  combinations  of  states  of  failed  functions 
X=triu( [ones(m,m) ; zeros (l,m)] ) ; 

7,  Create  a  matrix  S  where  each  row  represents  a  possible  combination  of  the 
7,  state  of  the  satellite . 
for  i=l : length(X) 

P=perms (X ( i , : ) ) ;  l  Create  permutations  of  having  i-1  function  down 

S=[S;  unique (P, ’rows’ )]  ;  7«  Delete  redundant  combinations  and  append  to  S 

end 

7,  Build  the  transistion  rate  generator  matrix. 

Q=[];  7o  Initialize  Q 

7,  Iterate  through  the  elements  of  S  finding  all  commutative  states.  When  a 
7,  state  is  found,  insert  the  corresponding  rate  into  the  Q  matrix, 
for  i=l:dimS 

for  j=l:dimS 

7,  Find  an  instance  where  the  number  of  functions  has  increased  or 
7,  decreased  by  1 . 
if  sum(S(i , : )==S(j , : ) )==(m-l) 
y=S(i, : ) — S ( j , : ) ; 

for  k=l:m  7*  Determine  which  function  changed, 

if  y (1 ,k)==l 

Q(i,j)=f(k);  7o  Function  failed, 

elseif  y(l,k)==-l 


A-17 


Q(i, j)=r(k) ; 


°/e  Function  was  repaired. 


end 

end 

end 

end 

end 

7.  Write  in  diagonals  on  condition  of  row  sums  equal  to  zero, 
for  i=l:dimS 

Q(i,i)=-sum(Q(i, : ))  ; 

end 


A. 2  Simulated  Availability 

function  [t,A_sim,  duration]  =  Sim_Examplel (t ,reps) 

7,  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR,  Penn  State  University 

7.  Date:  January  15,  2001 

7,  Revised  by:  Captain  Chris  Solo,  M.S.  candidate,  OR,  Air  Force  Institute  of  Technology 
7,  Date:  29  January  2004 

7,  Revised  by:  Captain  Cole  Gulyas,  M.S.  candidate,  OR,  Air  Force  Institute  of  Technology 
7,  Date:  10  November  2004 

7. 

7,  The  purpose  of  this  MATLAB  program  is  to  simulate  a  finite-state  semi-Markov  process. 

7,  The  process  is  simulated  in  order  to  validate  the  probability  mass  function  at  a  specific 

7,  time  gained  from  transient  analysis  a  satellite  system  modelled  as  an  SMP. 

7.  The  program  uses  function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 

7. 

7.  Input:  Q=Inf initesimal  generator  matrix,  t=time  associated  with  probability 

7,  mass  function,  a=initial  state  vector,  and  reps=number  of 

7.  simulation  repetitions 

7,  Output:  Measures  for  Example  1:  t=time,  A_sim=simulated  instantaneous  availability, 

7,  and  duration=processing  time. 

7ol-4  Block  A 
m=3 ; 
n=2~m; 

a=[l  0  0  0  0  0  0  0]  ; 
s=[l  0  0  0  0  0  0  0]  ; 
f =  [ . 0291  .0503  .0230]; 
r= [0  3.33  0] ; 
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7.  0/o5-16  Block  A 
7.  m=3; 

7.  n=2'~m; 

7.  a=[l  0  0  0  0  0  0  0]  ; 

'/,  s=[l  0  0  0  0  0  0  0]; 

7.  f  =  [  .0121  .0421  .0189]  ; 

7.  r=  [0  4.8216  0]  ; 

7.  7.17-20  Block  R 
7.  m=3; 

7.  n=2~m ; 

7.  a=  [ .  98  0  . 02  0  0  0  0  0]  ; 

7.  s=[l  0  0  0  0  0  0  0]  ; 

7.  f  =  [  .0121  .0303  .0143]  ; 

7.  r=  [0  4.9216  0]  ; 

7.  7.21-24  Block  R 
7.  m=3; 

7.  n=2~m ; 

7.  a=[l  0  0  0  0  0  0  0]  ; 

7.  s=[l  0  0  0  0  0  0  0]  ; 

7,  f  =  [ .  0103  .0282  .0118]  ; 

7.  r=  [0  4.9216  0]  ; 


Q=rate_matrix(m, f ,r) ; 


7.  Convert  mean  holding  times  for  each  state  from  lambda  form  to  mu 
K=-l . / diag(Q) ; 


tic 


7.  Start  timer 


7.  Probability  transition  matrix 
P=zeros (n,n) ; 
for  1=1 :n 

for  m=l:n 
if  l~=m 

P(l,m)=(Q(l,m)/ (— Q (1,1))) ; 

end 

end 

end 


statecounter  =  zeros (l,n);  7.  Counts  number  of  times  process  was  found  in  each  state 
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for  k  =  l:reps 


if  (mod(k,  10000)==0)  disp([’k  =  ’  num2str  (k)  ]  )  ;  end  7,  Display  number  of  reps 

Z  =  []; 

Z(l)  =  rando(a);  7*  Initial  state  of  the  environment  at  time  0 

newtime  =  0; 

7,  Specify  the  distribution  for  the  initial  state,  corresponding  to  vector  a 
totaltime  =  exprnd(K(Z(l) ) ) ;  %  ******  Time  spent  in  initial  state 

i=l ; 

while  (totaltime  <  t) 

Z(i+1)  =  rando(P(Z(i) , ;  l  Use  P  matrix  to  determine  next  state 

switch  Z(i+1) 
case  {1} 

newtime  =  exprnd(K(Z(i+l) ) ) ; 
case  {2} 

newtime  =  exprnd(K(Z(i+l) ) ) ; 
case  {3} 

newtime  =  exprnd(K(Z(i+l) ) ) ; 
case  {4} 

newtime  =  exprnd(K(Z(i+l) ) ) ; 
case  {5} 

newtime  =  exprnd(K(Z(i+l) ) ) ; 
case  {6} 

newtime  =  exprnd(K(Z(i+l) ) ) ; 
case  {7} 

newtime  =  exprnd(K(Z(i+l) ) ) ; 
case  {8} 

newtime  =  exprnd(K(Z(i+l) ) ) ; 

end 

totaltime  =  totaltime  +  newtime; 
i=i+l ; 

end 

statecounter (Z(end) )  =  statecounter (Z(end) )  +  1; 

end 

pmf=  statecounter . /reps ; 

A_sim=pmf *s ’ ; 
duration=toc 

function  [t,A_sim,  duration]  =  Sim_Example2 (t ,reps) 

7.  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR,  Penn  State  University 
7,  Date:  January  15,  2001 
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7,  Revised  by:  Captain  Chris  Solo,  M.S.  candidate,  OR,  Air  Force  Institute  of  Technology 
7*  Date:  29  January  2004 

7,  Revised  by:  Captain  Cole  Gulyas,  M.S.  candidate,  OR,  Air  Force  Institute  of  Technology 
7.  Date:  6  January  2005 

7. 

7,  The  purpose  of  this  MATLAB  program  is  to  simulate  a  finite-state  semi-Markov  process. 

7,  The  process  is  simulated  in  order  to  demonstrate  the  probability  mass  function  at  a  specific 
7,  time  gained  from  transient  analysis  of  a  satellite  system  modelled  as  an  SMP. 

7,  The  program  uses  function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 

7. 

7.  Input:  t=time  associated  with  probability  mass  function 

7,  reps=number  of  simulation  repetitions 

7,  Output:  Measures  for  Example  2:  t=time,  A_sim=simulated  instantaneous  availability, 

7,  and  duration=processing  time. 


m=4 ; 
n=2~m; 

a=  [1  00000000000000  0]; 
s=  [1  01110000000000  0]; 


K=  [3.8 

.37 

0 

7.  Gamma 

2.1 

.50 

0 

7.  Gamma 

2.2 

.41 

0 

7,  Gamma 

3.2 

.42 

0 

7.  Gamma 

1.2 

.22 

0 

7.  Weibull 

1.8 

.42 

0 

7.  Weibull 

1.9 

.31 

0 

7.  Weibull 

1.2 

.33 

0 

7.  Weibull 

1.8 

2.1 

2.3 

7,  Triangular 

1.61  2.2 

3.2 

7.  Triangular 

1.1 

1.83 

2.4 

7,  Triangular 

1.7 

2.5 

4.51 

7.  Triangular 

.05 

2.3 

0 

7,  Uniform 

.  13 

2.1 

0 

7,  Uniform 

.03 

2.4 

0 

7,  Uniform 

.51 

.93 

0]; 

7.  Uniform 

7,  Mean  holding  times 

for  each  state 

mu=  []  ; 
for  i=l:n 

mu=  [mu;  (1/K(i ,2) ) *gamma( (K(i, 1)+1)/K(i , 1) )] ;  7»k/lamda 

end 


R=[-l/mu(l)  . 2/mu(l)  .3/mu(l)  .3/mu(l)  .2/mu(l)  00000000000 
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. 2/mu(2)  -l/mu(2)  000  .4/mu(2)  .2/mu(2)  .2/mu(2)  00000000 
0  0  -l/mu(3)  0  0  . 2/mu (3)  0  0  .4/mu(3)  .4/mu(3)  000000 
000  -l/mu(4)  0  0  . 5/mu(4)  0  .3/mu(4)  0  .2/mu(4)  00000 
0000  -l/mu(5)  0  0  .2/mu(5)  0  .7/mu(5)  .l/mu(5)  00000 
0  0  . 4/mu(6)  0  0  -l/mu(6)  00000  .3/mu(6)  .3/mu(6)  000 
000  . 5/mu(7)  0  0  -l/mu(7)  0000  .2/mu(7)  0  .3/mu(7)  0  0 
0000  . 4/mu(8)  0  0  -l/mu(8)  0000  .5/mu(8)  .l/rau(8)  0  0 
00000000  -l/mu(9)  0  0  .7/mu(9)  0  0  .3/mu(9)  0 
000000000  -l/mu(10)  0  0  .8/mu(10)  0  .2/mu(10)  0 
0000000000  -l/mu(ll)  0  0  .6/mu(ll)  .4/mu(ll)  0 
00000000  . 2/mu(12)  0  0  -l/mu(12)  000  .8/mu(12) 
000000000  . 3/mu(13)  0  0  -l/mu(13)  0  0  .7/mu(13) 
0000000000  . 3/mu(14)  0  0  -l/mu(14)  0  .7/mu(14) 
00000000000000  -l/mu(15)  l/mu(15) 
000000000000000  eps] ; 


7.  Convert  Gamma  parameters  from  f  (x I k, lambda)  to  f(x|k,beta) 
7.  in  order  to  use  the  matlab  function  weibrnd 
for  i=l:4 

K(i,:)=[K(i,l)  1/K(i ,2)  0]; 


end 


7.  Convert  Weibull  parameters  from  f  (x  I  alpha,  lambda)  to  f(x|a,b) 
7.  in  order  to  use  the  matlab  function  weibrnd 
for  i=5:8 

K(i,:)=[K(i,2)~K(i,l)  K(i,l)  0]; 


end 


tic  7.  Start  timer 

7,  Probability  transition  matrix 
P=zeros (n,n) ; 
for  1=1 :n 

for  m=l:n 
if  l~=m 

P(l,m)=(R(l,m)/(-R(l,l))); 

end 

end 

end 


statecounter  =  zeros (l,n);  7«  Counts  number  of  times  process  was  found  in  each  state 


for  k  =  l:reps 
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if  (mod(k,  1000)==0)  disp([’k  =  ’  num2str (k) ]  )  ;  end  "l  Display  number  of  reps 

Z  =  []; 

Z(l)  =  rando(a);  %  Initial  state  of  the  environment  at  time  0 

newtime  =  0; 

°/«  Specify  the  distribution  for  the  initial  state,  corresponding  to  vector  a 
totaltime  =  gamrnd(K(Z(l)  ,  1)  ,K(Z(1)  ,2) )  ;  °/0  ******  Time  spent  in  initial  state 

i=l ; 

while  (totaltime  <  t) 

Z(i+1)  =  rando(P(Z(i)  ,  : ) )  ;  °/«  Use  P  matrix  to  determine  next  state 

switch  Z(i+1) 
case  {1} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {2} 

newtime  =  gamrnd(K(Z(i+l) , 1) , l/K(Z(i+l) ,2) ) ; 
case  {3} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {4} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {5} 

newtime  =  weibrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {6} 

newtime  =  weibrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {7} 

newtime  =  weibrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {8} 

newtime  =  weibrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {9} 

newtime  =  trirnd(K(Z(i+l) ,1) ,K(Z(i+l) ,2) ,K(Z(i+l) ,3)) ; 
case  {10} 

newtime  =  trirnd(K(Z(i+l) ,1) ,K(Z(i+l) ,2) ,K(Z(i+l) ,3)) ; 
case  {11} 

newtime  =  trirnd(K(Z(i+l) ,1) ,K(Z(i+l) ,2) ,K(Z(i+l) ,3)) ; 
case  {12} 

newtime  =  trirnd(K(Z(i+l) ,1) ,K(Z(i+l) ,2) ,K(Z(i+l) ,3)) ; 
case  {13} 

newtime  =  unif rnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {14} 

newtime  =  unif rnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {15} 

newtime  =  unif rnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {16} 
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newtime  =  unif rnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 


end 

totaltime  =  totaltime  +  newtime; 
i=i+l ; 

end 

statecounter (Z(end) )  =  statecounter (Z(end) )  +  1; 

end 

pmf=  statecounter . /reps ; 

A_sim=pmf *s ’ ; 
duration=toc 


function  [t,A_sim,  duration]  =  Sim_Example3_l (t ,reps) 

7.  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR,  Penn  State  University 

7*  Date:  January  15,  2001 

7,  Revised  by:  Captain  Chris  Solo,  M.S.  candidate,  OR,  Air  Force  Institute  of  Technology 
7,  Date:  29  January  2004 

7,  Revised  by:  Captain  Cole  Gulyas,  M.S.  candidate,  OR,  Air  Force  Institute  of  Technology 
7,  Date:  6  January  2005 

7. 


7.  The  purpose  of  this  MATLAB  program  is  to  simulate  a  finite-state  semi-Markov  process. 

7,  The  process  is  simulated  in  order  to  demonstrate  the  probability  mass  function  at  a  specific 
7,  time  gained  from  transient  analysis  of  a  satellite  system  modelled  as  an  SMP. 

7,  The  program  uses  function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 


7. 

7,  Input : 

7. 

7,  Output : 
7. 


t=time  associated  with  probability  mass  function 
reps=number  of  simulation  repetitions 

Measures  for  Example  3,  satellite  1:  t=time,  A_sim=simulated  instantaneous 
availability,  and  duration=processing  time. 


m=4 ; 
n=2~m ; 

a=  [1  00000000000000  0]; 
s=  [1  00000000000000  0]; 


[5.8 

.37  0 

7,  Gamma 

4.1 

6.80  0 

7.  Gamma 

1.6 

2.31  0 

7.  Gamma 

1.2 

1.13  0 

7,  Gamma 

2  7. 

.41  0 

7,  Erlang 

5  6. 

.42  0 

7,  Erlang 

1.2 

1.22  0 

7.  Weibull 

1.8 

2.2  0 

7.  Weibull 

.05 

1.4  0 

7.  Uniform 
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.13  1.1  0 

7.  Uniform 

.02  .5  1.3 

7.  Triangular 

.13  .4  1.2 

7,  Triangular 

.09  .3  1.1 

7,  Triangular 

11.6  0 

7«  Exp 

11.4  0 

7.  Exp 

1  1.8  0  ] ; 

7«  Exp 

7.  Mean  holding  times  for  each  state 
mu=  []  ; 
for  i=l:6 

mu=  [mu;  K(i ,  1) /K(i ,2)]  ;  70Gamma  &  Erlang 

end 

for  i=7:8 

mu=  [mu;  (1/K(i ,2) ) *gamma( (K(i , 1)+1)/K(i , 1) )] ;  7,Weibull 

end 

for  i=9:10 

mu=  [mu;  (K(i ,  1)  +K(i , 2) ) / 2]  ;  70Uniform 

end 


for  i=ll:13 

mu=  [mu;  (K(i,l)+K(i,2)+K(i,3))/3]  ;  70Triangular 

end 

for  i=14:16 

mu=  [mu;  K(i ,  1) /K(i ,2)]  ;  70Exp 

end 


7*  Rate  matrix  built  with  all  types  of  repairs  possible 
R=[-l/mu(l)  . 2/mu(l)  .3/mu(l)  .3/mu(l)  .2/mu(l)  00000000000 
. 4/mu(2)  -l/mu(2)  000  .2/mu(2)  .2/mu(2)  .2/mu(2)  00000000 
. 5/mu (3)  0  -l/mu(3)  0  0  .l/mu(3)  0  0  .2/mu(3)  .2/mu(3)  000000 
. 2/mu (4)  0  0  -l/mu(4)  0  0  .4/mu(4)  0  .2/mu(4)  0  .2/mu(4)  00000 
.4/mu(5)  000  -l/mu(5)  0  0  .l/mu(5)  0  .4/mu(5)  .l/mu(5)  00000 
0  . 3/mu(6)  . 3/mu(6)  0  0  -l/mu(6)  00000  .l/mu(6)  .3/mu(6)  000 
0  . 2/mu(7)  0  . 4/mu(7)  0  0  -l/mu(7)  0000  .2/mu(7)  0  .2/mu(7)  0  0 
0  . 3/mu(8)  0  0  . 3/mu(8)  0  0  -l/mu(8)  0000  .3/mu(8)  .l/mu(8)  0  0 
0  0  . 3/mu(9)  . 3/mu(9)  0000  -l/mu(9)  0  0  .3/mu(9)  0  0  .l/mu(9)  0 
0  0  . 2/mu(10)  0  . 3/mu(10)  0000  -l/mu(10)  0  0  .3/mu(10)  0  .2/mu(10)  0 

000  . 3/mu(ll)  . l/mu(ll)  00000  -l/mu(ll)  0  0  .3/mu(ll)  .3/mu(ll)  0 

00000  . 2/mu(12)  .3/mu(12)  0  .2/mu(12)  0  0  -l/mu(12)  000  .3/mu(12) 
00000  . 3/mu (13)  0  .3/mu(13)  0  .l/mu(13)  0  0  -l/mu(13)  0  0  .3/mu(13) 

000000  . 2/mu(14)  .4/mu(14)  0  0  .3/mu(14)  0  0  -l/mu(14)  0  .l/mu(14) 

00000000  . 2/mu(15)  .3/mu(15)  .3/mu(15)  000  -l/mu(15)  .2/mu(15) 
00000000000  . 4/mu(16)  .4/mu(16)  .l/mu(16)  .l/mu(16)  -l/mu(16)] ; 
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7,  Convert  Gamma  parameters  from  f  (x I k, lambda)  to  f(x|k,beta) 

7,  in  order  to  use  the  matlab  function  gamrnd 
for  i=l:6 

K(i,:)=[K(i,l)  1/K(i ,2)  0]; 

end 

7,  Convert  Weibull  parameters  from  f  (x  I  alpha, lambda)  to  f(x|a,b) 

7,  in  order  to  use  the  matlab  function  weibrnd 
for  i=7:8 

K(i,:)=[K(i,2)~K(i,l)  K(i,l)  0]; 

end 

7.  Convert  Exp  parameters  from  f(x  I  lambda)  to  f(x|beta) 

7.  in  order  to  use  the  matlab  function  gamrnd 
for  i=14:16 

K(i, :)=[K(i,l)  1/K(i,2)  0]; 

end 

tic  7.  Start  timer 

7,  Probability  transition  matrix 
P=zeros (n,n) ; 
for  1=1 :n 

for  m=l:n 
if  l~=m 

P(l,m)=(R(l,m)/(-R(l , 1) ) ) ; 

end 

end 

end 

statecounter  =  zeros (l,n);  7«  Counts  number  of  times  process  was  found  in  each  state 

for  k  =  l:reps 

if  (mod(k, 1000)==0)  disp([’k 
Z  =  []; 

Z(l)  =  rando(a) ; 
newtime  =  0; 

7.  Specify  the  distribution  for  the  initial  state,  corresponding  to  vector  a 
totaltime  =  gamrnd(K(Z(l) , 1) ,K(Z(1) ,2) ) ;  %  ******  Time  spent  in  initial  state 

i=l ; 


=  ’  num2str  (k)  ]  )  ;  end  7«  Display  number  of  reps 


7,  Initial  state  of  the  environment  at  time  0 
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while  (totaltime  <  t) 

Z(i+1)  =  rando(P(Z(i)  ,  : ) )  ;  °/«  Use  P  matrix  to  determine  next  state 

switch  Z(i+1) 
case  {1} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {2} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {3} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {4} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) , 2) ) ; 
case  {5} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {6} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {7} 

newtime  =  weibrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {8} 

newtime  =  weibrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {9} 

newtime  =  unif rnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {10} 

newtime  =  unif rnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {11} 

newtime  =  trirnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ,K(Z(i+l) ,3) ) ; 
case  {12} 

newtime  =  trirnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ,K(Z(i+l) ,3) ) ; 
case  {13} 

newtime  =  trirnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ,K(Z(i+l) ,3) ) ; 
case  {14} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {15} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {16} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 

end 

totaltime  =  totaltime  +  newtime; 
i=i+l ; 

end 

statecounter (Z(end) )  =  statecounter (Z(end) )  +  1; 

end 

pmf=  statecounter . /reps ; 
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A_sim=pmf *s ’ ; 
duration=toc 

function  [t,A_sim,  duration]  =  Sim_Example3_2 (t ,reps) 

7,  Orig  Author:  Jeffrey  P.  Kharoufeh,  Ph.D.  candidate,  IE  &  OR,  Penn  State  University 
7*  Date:  January  15,  2001 

7,  Revised  by:  Captain  Chris  Solo,  M.S.  candidate,  OR,  Air  Force  Institute  of  Technology 
7,  Date:  29  January  2004 

7,  Revised  by:  Captain  Cole  Gulyas,  M.S.  candidate,  OR,  Air  Force  Institute  of  Technology 

7,  Date:  6  January  2005 

7. 

7,  The  purpose  of  this  MATLAB  program  is  to  simulate  a  finite-state  semi-Markov  process. 

7.  The  process  is  simulated  in  order  to  demonstrate  the  probability  mass  function  at  a  specific 

7,  time  gained  from  transient  analysis  of  a  satellite  system  modelled  as  an  SMP. 

7,  The  program  uses  function  "rando"  in  order  to  select  the  next  state  after  a  state  transition. 

7. 

7,  Input:  t=time  associated  with  probability  mass  function 

7,  reps=number  of  simulation  repetitions 

7,  Output:  Measures  for  Example  3,  satellite  2:  t=time,  A_sim=simulated  instantaneous 

7.  availability,  and  duration=processing  time. 

m=3 ; 
n=2~m; 

a=[l  0  0  0  0  0  0  0]  ; 
s=[l  0  0  0  0  0  0  0]  ; 


[8.7  1.07  0 

7o  Gamma 

3.1  5.80  0 

7.  Gamma 

2  6.41  0 

7.  Erlang 

5  6 . 42  0 

7.  Erlang 

1.2  1.22  0 

7.  Weibull 

1.8  2.2  0 

7.  Weibull 

.02  .5  1.3 

7.  Triangular 

1  1.8  0  ] ; 

7*  Exp 

7,  Mean  holding  times  for  each  state 
mu=  []  ; 
for  i=l:4 

mu= [mu;  K(i, 1)/K(i ,2)] ; 

end 

for  i=5:6 

mu= [mu;  (1/K(i ,2) ) *gamma( (K(i , 1)+1) /K(i , 1) )] ; 

end 

for  i=7 


70Gamma  &  Erlang 


7oWeibull 
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mu=[mu;  (K(i , 1)+K(i , 2)+K(i , 3) ) / 3] ; 

end 

for  i=8 

mu= [mu;  K(i , 1) /K(i ,2)] ; 

end 

7,  Rate  matrix  built  with  all  types  of  repairs  possible 
7*  R  matrix  (observed  rates  of  state  transitions) 

R= [  -l/mu(l)  . 5/mu(l)  .3/mu(l)  .2/mu(l)  0000 
.4/mu(2)  -l/mu(2)  0  0  .3/mu(2)  .3/mu(2)  0  0 
.3/mu (3)  0  -1/mu (3)  0  .3/mu (3)  0  .4/mu (3)  0 
.3/mu(4)  0  0  -l/mu(4)  0  .5/mu(4)  .2/mu(4)  0 
0  .5/mu(5)  .4/mu(5)  0  -l/mu(5)  0  0  .l/mu(5) 

0  .4/mu(6)  0  .3/mu(6)  0  -l/mu(6)  0  .3/mu(6) 

0  0  .4/mu(7)  . 3/mu (7)  0  0  -l/mu(7)  .3/mu(7) 

0000  .3/mu(8)  .3/mu(8)  .4/mu(8)  -l/mu(8)] ; 

7.  Convert  Gamma  parameters  from  f (x I k, lambda)  to  f(x|k,beta) 

7.  in  order  to  use  the  matlab  function  weibrnd 
for  i=l:4 

K(i,:)=[K(i,l)  1/K(i , 2)  0]; 

end 

7,  Convert  Weibull  parameters  from  f  (x  I  alpha,  lambda)  to  f(x|a,b) 
7,  in  order  to  use  the  matlab  function  weibrnd 
for  i=5:6 

K(i,:)=[K(i,2)~K(i,l)  K(i,l)  0]; 

end 

7,  Convert  Gamma  parameters  from  f  (x I k, lambda)  to  f(x|k,beta) 

7.  in  order  to  use  the  matlab  function  weibrnd 
for  i=8 

K(i , : )  =  [K(i , 1)  1/K(i,2)  0]; 

end 

tic  7.  Start  timer 

7,  Probability  transition  matrix 
P=zeros (n,n) ; 
for  1=1 :n 

for  m=l:n 
if  l~=m 

P(1 ,m)=(R(l,m)/(-R(l , 1) ) ) ; 


7»Tri  angular 


7oExp 
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end 


end 


end 

statecounter  =  zeros (l,n);  7«  Counts  number  of  times  process  was  found  in  each  state 

for  k  =  l:reps 

if  (mod(k, 1000)==0)  disp([’k 
Z  =  []; 

Z(l)  =  rando(a) ; 
newtime  =  0; 

7.  Specify  the  distribution  for  the  initial  state,  corresponding  to  vector  a 
totaltime  =  gamrnd(K(Z(l)  ,  1)  ,K(Z(1)  ,2) )  ;  °/0  ******  Time  spent  in  initial  state 

i=l ; 

while  (totaltime  <  t) 

Z(i+1)  =  rando(P(Z(i) ,  ;  °/«  Use  P  matrix  to  determine  next  state 

switch  Z(i+1) 
case  {1} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) , 2) ) ; 
case  {2} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {3} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {4} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {5} 

newtime  =  weibrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {6} 

newtime  =  weibrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 
case  {7} 

newtime  =  trirnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ,K(Z(i+l) ,3) ) ; 
case  {8} 

newtime  =  gamrnd(K(Z(i+l) , 1) ,K(Z(i+l) ,2) ) ; 

end 

totaltime  =  totaltime  +  newtime; 
i=i+l ; 

end 

statecounter (Z (end) )  =  statecounter (Z (end) )  +  1; 

end 

pmf=  statecounter . /reps ; 


=  3  num2str (k)]  )  ;  end  °/0  Display  number  of  reps 


7*  Initial  state  of  the  environment  at  time  0 
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A_sim=pmf *s ’ ; 
duration=toc 

function  [index]  =  rando(p) 

7,  rando.m  generates  a  discrete  random  variable  in  S={1 ,2, .  .  .  ,n}  given 
7«  a  distribution  vector  p  =  [pi  p2  ...  pn]  . 

u  =  rand; 
i  =  1; 
s  =  p(l) ; 

while  ((u  >  s)  &  (i  <  length(p))), 
i=i+l; 
s=s+p(i) ; 
end 

index=i ; 

function  n  =  trirnd(a,m,b) 

7*  Author:  Capt  Cole  Gulyas 

7,  Date:  15  Nov  04 

7.  Input:  a=min,  m=mode,  b=max 

°/«  Output:  Random  number  from  triangular  (a,  m,b) 

format  long 

d=(m-a)/(b-a) ; 
r=rand(l) ; 
if  r<d 

n=sqrt (r* (b-a) * (m-a) )  +  a; 

else 

n=b  -  sqrt ( (1-r) * (b-a) * (b-m) ) ; 

end 
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